From 36af3e3f73154e5e14c4745304bd409eaf265899 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 10 Jun 2026 15:07:30 +0200 Subject: [PATCH 1/3] Clean up image simulation code. --- app/pywmtk/pywmtk.cpp | 15 - cmake/recipes/libigl.cmake | 2 +- .../image_simulation/CMakeLists.txt | 4 - .../image_simulation/EmbedSurface.cpp | 456 +----------- .../image_simulation/EmbedSurface.hpp | 45 -- .../image_simulation/ImageSimulationMesh.cpp | 48 -- .../image_simulation/ImageSimulationMesh.h | 29 - .../ImageSimulationMeshTri.cpp | 85 +-- .../ImageSimulationMeshTri.hpp | 1 - .../components/image_simulation/Parameters.h | 49 +- .../components/image_simulation/Smooth.cpp | 34 +- .../VolumemesherInsertion.cpp | 703 ------------------ .../image_simulation/extract_soup.cpp | 395 ---------- .../image_simulation/extract_soup.hpp | 34 - .../image_simulation/image_simulation.cpp | 381 +++------- .../image_simulation_spec.hpp | 38 - .../image_simulation_spec.json | 38 - .../image_simulation/read_image_msh.cpp | 87 +-- .../image_simulation/read_image_msh.hpp | 10 - 19 files changed, 158 insertions(+), 2296 deletions(-) delete mode 100644 components/image_simulation/wmtk/components/image_simulation/extract_soup.cpp delete mode 100644 components/image_simulation/wmtk/components/image_simulation/extract_soup.hpp diff --git a/app/pywmtk/pywmtk.cpp b/app/pywmtk/pywmtk.cpp index 335b137eb1..20693c29af 100644 --- a/app/pywmtk/pywmtk.cpp +++ b/app/pywmtk/pywmtk.cpp @@ -11,10 +11,6 @@ // components #include "../components_include.hpp" -// std::map> components_map; -//// include auto-generated map -//#include "../components_map.hpp" - namespace py = pybind11; using namespace pybind11::literals; @@ -47,16 +43,5 @@ void wmtk_wrapper(const py::dict& obj) PYBIND11_MODULE(pywmtk, m, py::mod_gil_not_used()) { m.doc() = "Python bindings for the Wildmeshing-Toolkit"; // optional module docstring - // m.def("get_metrics", &meme::get_metrics, "Get mesh metrics"); - // m.def("get_metric_names", &meme::get_metrics_names, "Get names for all mesh metrics"); - // m.def("get_metrics_per_tri", &meme::get_metrics_per_tri, "Get mesh metrics per triangle"); - // m.def( - // "get_metric_names_per_tri", - // &meme::get_metrics_names_per_tri, - // "Get names for all per-triangle metrics"); - // m.def( - // "get_relative_edge_lenghts", - // &meme::get_relative_edge_lengths, - // "Get edge lengths realtive to bbox diagonal"); m.def("wmtk", &wmtk_wrapper, "Wildmeshing-Toolkit application"); } \ No newline at end of file diff --git a/cmake/recipes/libigl.cmake b/cmake/recipes/libigl.cmake index 916903c10f..be4f69d522 100644 --- a/cmake/recipes/libigl.cmake +++ b/cmake/recipes/libigl.cmake @@ -12,7 +12,7 @@ CPMAddPackage( GIT_TAG 3ea7f9480967fcf6bf02ce9b993c0ea6d2fc45f6 OPTIONS LIBIGL_PREDICATES ON - LIBIGL_COPYLEFT_TETGEN ON + # LIBIGL_COPYLEFT_TETGEN ON ) # include(eigen) FetchContent_MakeAvailable(libigl) \ No newline at end of file diff --git a/components/image_simulation/wmtk/components/image_simulation/CMakeLists.txt b/components/image_simulation/wmtk/components/image_simulation/CMakeLists.txt index a0e27f3b80..6a995db2f5 100644 --- a/components/image_simulation/wmtk/components/image_simulation/CMakeLists.txt +++ b/components/image_simulation/wmtk/components/image_simulation/CMakeLists.txt @@ -13,9 +13,6 @@ set(SRC_FILES ImageSimulationMesh.cpp ImageSimulationMesh.h - extract_soup.hpp - extract_soup.cpp - EmbedSurface.hpp EmbedSurface.cpp @@ -54,7 +51,6 @@ target_link_libraries(wmtk_${COMPONENT_NAME} PUBLIC wmtk::shortest_edge_collapse VolumeRemesher::VolumeRemesher jse::jse - igl_copyleft::tetgen ) ###################### diff --git a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp index 3cb1e2bba0..2ddce88bc6 100644 --- a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp @@ -5,7 +5,6 @@ #include #include #include -#include #include #include #include @@ -28,8 +27,6 @@ #include -#include "extract_soup.hpp" - namespace { std::vector> triangulate_polygon_face(std::vector points) { @@ -183,18 +180,6 @@ void delaunay_box_mesh( points.size()); } -void embed_surface( - const MatrixXd& V_surface, - const MatrixXi& F_surface, - const MatrixXd& V_vol, - const MatrixXi& T_vol, - MatrixXr& V_emb, - MatrixXi& T_emb) -{ - MatrixXi F_on_surface; - embed_surface(V_surface, F_surface, V_vol, T_vol, V_emb, T_emb, F_on_surface); -} - void embed_surface( const MatrixXd& V_surface, const MatrixXi& F_surface, @@ -564,297 +549,6 @@ void embed_surface( std::cout << "v on surface: " << v_on_input.size() << std::endl; } -void tag_tets_from_image( - const std::string& filename, - const Matrix4d& xyz2ijk, - const MatrixXd& V, - const MatrixXi& T, - VectorXi& T_tags) - -{ - logger().info("Tag tets"); - T_tags.resize(T.rows(), 1); - T_tags.setZero(); - - std::vector>> volumetric_data; - read_array_data_ascii(volumetric_data, filename); - - tag_tets_from_image(volumetric_data, xyz2ijk, V, T, T_tags); -} - -void tag_tets_from_image( - const ImageData& data, - const Matrix4d& xyz2ijk, - const MatrixXd& V, - const MatrixXi& T, - VectorXi& T_tags) -{ - T_tags.resize(T.rows(), 1); - T_tags.setZero(); - - for (size_t i = 0; i < T.rows(); ++i) { - const Vector3d v0 = V.row(T(i, 0)); - const Vector3d v1 = V.row(T(i, 1)); - const Vector3d v2 = V.row(T(i, 2)); - const Vector3d v3 = V.row(T(i, 3)); - - Vector3d center = (v0 + v1 + v2 + v3) * 0.25; - // map from geometric position to index space, i.e., undo dimension multiplication - //center[0] /= dimensions[0]; - //center[1] /= dimensions[1]; - //center[2] /= dimensions[2]; - { - Vector4d ch = to_homogenuous(center); - ch = xyz2ijk * ch; - center = from_homogenuous(ch); - } - const int idx_0 = std::floor(center.x()); - const int idx_1 = std::floor(center.y()); - const int idx_2 = std::floor(center.z()); - // mesh.m_tet_attribute[t.tid(mesh)].tag = idx_0; - if (idx_0 >= 0 && idx_0 < data.size() && idx_1 >= 0 && idx_1 < data[0].size() && - idx_2 >= 0 && idx_2 < data[0][0].size()) { - // for tag - int64_t intValue = data[idx_0][idx_1][idx_2]; - T_tags[i] = intValue; - } - } -} - -void tag_tets_from_images( - const std::vector& data, - const Matrix4d& xyz2ijk, - const MatrixXd& V, - const MatrixXi& T, - MatrixSi& T_tags) -{ - T_tags.resize(T.rows(), data.size()); - T_tags.setZero(); - - for (size_t i = 0; i < data.size(); ++i) { - VectorXi t; - tag_tets_from_image(data[i], xyz2ijk, V, T, t); - for (size_t j = 0; j < t.size(); ++j) { - if (t[j] == 0) { - continue; - } - T_tags.coeffRef(i, j) = t[j]; - } - } - T_tags.makeCompressed(); -} - -EmbedSurface::EmbedSurface(const std::vector& img_filenames, const Matrix4d& ijk2xyz) - : m_img_filenames(img_filenames) - , m_ijk2xyz(ijk2xyz) - , m_xyz2ijk(ijk2xyz.inverse()) -{ - m_img_datas.resize(m_img_filenames.size()); - for (size_t i = 0; i < m_img_filenames.size(); ++i) { - if (!std::filesystem::exists(m_img_filenames[i])) { - log_and_throw_error("Image {} does not exist", m_img_filenames[i]); - } - read_array_data_ascii(m_img_datas[i], m_img_filenames[i]); - MatrixXd V_single; - MatrixXi F_single; - MatrixXi F_tags_single; - extract_triangle_soup_from_image(m_img_datas[i], F_single, V_single, F_tags_single); - - assert(V_single.cols() == 3); - assert(F_single.cols() == 3); - - const size_t nV_old = m_V_surface.rows(); - const size_t nF_old = m_F_surface.rows(); - - m_V_surface.conservativeResize(m_V_surface.rows() + V_single.rows(), 3); - m_V_surface.block(nV_old, 0, V_single.rows(), 3) = V_single; - - F_single.array() += nV_old; - m_F_surface.conservativeResize(m_F_surface.rows() + F_single.rows(), 3); - m_F_surface.block(nF_old, 0, F_single.rows(), 3) = F_single; - - m_F_tags_surface.reserve(m_F_tags_surface.size() + F_tags_single.rows()); - for (size_t j = 0; j < F_tags_single.rows(); ++j) { - std::vector> tags = { - std::make_pair(i, (size_t)F_tags_single(j, 0)), - std::make_pair(i, (size_t)F_tags_single(j, 1))}; - m_F_tags_surface.push_back(tags); - } - } - - - // process triangle soup - - MatrixXd V; - MatrixXi F; - VectorXi _I; - - // remove unreferenced vertices - igl::remove_unreferenced(m_V_surface, m_F_surface, V, F, _I); - - if (V.rows() == 0 || F.rows() == 0) { - log_and_throw_error("Empty Input"); - } - - // const std::pair box_minmax = - // std::pair(V.colwise().minCoeff(), V.colwise().maxCoeff()); - // double diag = (box_minmax.first - box_minmax.second).norm(); - - // using the same error tolerance as in tetwild - Eigen::VectorXi SVI, SVJ, SVK; - Eigen::MatrixXd temp_V = V; // for STL file - - const Vector4d single_voxel_max = ijk2xyz * Vector4d::Ones(); - const Vector4d single_voxel_min = ijk2xyz * Vector4d(0, 0, 0, 1); - //// ??? this might be wrong. I think we are still in voxel space here, so eps should be a fixed - //// value like 0.1. - // double eps = (from_homogenuous(single_voxel_max) - from_homogenuous(single_voxel_min)) - // .cwiseAbs() - // .minCoeff() * - // 0.01; - // if (eps <= 0) { - // logger().warn("EPS = {}, ijk_to_ras matix might be broken! Changing eps to 1e-4", eps); - // eps = 1e-4; - // } - double eps = 0.1; - - igl::remove_duplicate_vertices(temp_V, eps, V, SVI, SVJ); - for (int i = 0; i < F.rows(); i++) { - for (int j : {0, 1, 2}) { - F(i, j) = SVJ[F(i, j)]; - } - } - - // remove duplicate faces and merge tag data - { - MatrixXi F1 = F; - VectorXi J; - resolve_duplicated_faces(F1, F, J); - - std::vector>> tags_surface; - tags_surface.resize(F.rows()); - for (size_t i = 0; i < J.rows(); ++i) { - if (i == J[i]) { - tags_surface[i] = m_F_tags_surface[i]; - } else { - tags_surface[J[i]].insert( - tags_surface[J[i]].end(), - m_F_tags_surface[i].begin(), - m_F_tags_surface[i].end()); - } - } - m_F_tags_surface = tags_surface; - } - wmtk::logger().info("after remove duplicate v#: {} f#: {}", V.rows(), F.rows()); - - // uniform laplacian smoothing - if (m_smooth_surface) { - logger().info("Taubin smoothing on input"); - // igl::writeOFF("debug_in.off", V, F); - const MatrixXd V_orig = V; - MatrixXd V_buf = V; - - const double lambda = 0.5; - const double my = -0.53; - // const double clamp_to = 0.249; - const double clamp_to = std::sqrt(0.1249); - - std::vector> A; - igl::adjacency_list(F, A); - - std::vector vertex_is_nonmanifold(V.rows(), false); - { - // find vertices on non-manifold edges - std::map edge_face_counter; - for (int i = 0; i < F.rows(); ++i) { - for (int j = 0; j < 3; ++j) { - const simplex::Edge e(F(i, j), F(i, (j + 1) % 3)); - auto [it, suc] = edge_face_counter.try_emplace(e, 0); - it->second++; - } - } - for (const auto& [e, fs] : edge_face_counter) { - if (fs != 2) { - vertex_is_nonmanifold[e.vertices()[0]] = true; - vertex_is_nonmanifold[e.vertices()[1]] = true; - } - } - } - - for (int i = 0; i < V.rows(); ++i) { - const auto& neighs = A[i]; - if (!vertex_is_nonmanifold[i]) { - continue; - } - std::vector n; - n.reserve(neighs.size()); - for (const int neigh : neighs) { - if (vertex_is_nonmanifold[neigh]) { - n.push_back(neigh); - } - } - A[i] = n; - } - - auto smooth_step = [&](const double factor) { - for (int i = 0; i < V.rows(); ++i) { - const auto& neighs = A[i]; - if (neighs.empty()) { - continue; - } - - Vector3d Lp = Vector3d::Zero(); - for (const int j : neighs) { - Lp += V.row(j); - } - Lp /= neighs.size(); - Lp -= V.row(i); - - V_buf.row(i) = V.row(i).transpose() + factor * Lp; - } - std::swap(V, V_buf); - }; - - auto clamp_step = [&]() { - for (int i = 0; i < V.rows(); ++i) { - const Vector3d& p0 = V_orig.row(i); - const Vector3d& p1 = V.row(i); - Vector3d dir = (p1 - p0); - // dir = dir.cwiseMax(-clamp_to).cwiseMin(clamp_to); - if (dir.norm() > clamp_to) { - dir = dir.normalized() * clamp_to; - } - V.row(i) = p0 + dir; - } - }; - - for (size_t n = 0; n < 100; ++n) { - smooth_step(lambda); - smooth_step(my); - clamp_step(); - } - // igl::writeOFF("debug_out.off", V, F); - } - - std::vector verts; - std::vector> tris; - verts.resize(V.rows()); - tris.resize(F.rows()); - wmtk::eigen_to_wmtk_input(verts, tris, V, F); - - // apply dimensions - logger().info("Convert from image to world coordinates (ijk to xyz)"); - for (Vector3d& v : verts) { - Vector4d vh = to_homogenuous(v); - vh = m_ijk2xyz * vh; - v = from_homogenuous(vh); - } - logger().info("Finished conversion."); - - V_surf_from_vector(verts); - F_surf_from_vector(tris); -} - EmbedSurface::EmbedSurface( const std::vector& img_filenames, const std::vector& img_transform, @@ -1034,12 +728,10 @@ bool EmbedSurface::embed_surface(const bool flood_fill) } // add tags - if (!Fs.empty()) { - tag_from_winding_number(); - } else { - assert(!m_img_datas.empty()); - tag_tets_from_images(m_img_datas, m_xyz2ijk, m_V_emb, m_T_emb, m_T_tags); + if (Fs.empty()) { + log_and_throw_error("No input surface to embed"); } + tag_from_winding_number(); /** * Cluster tags by flood-filling regions that are bounded by the surface. All tags within one @@ -1131,148 +823,6 @@ bool EmbedSurface::embed_surface(const bool flood_fill) return all_rounded; } -bool EmbedSurface::embed_surface_tetgen() -{ - logger().info("Embed with TetGen"); - - std::shared_ptr ptr_env; - { - const auto v_simplified = V_surf_to_vector(); - - std::vector 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(); - ptr_env->use_exact = true; - 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 points(8); - for (int i = 0; i < 8; i++) { - Vector3d& p = points[i]; - std::bitset 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 N; // number of grid points per dimension - std::array 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, 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 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; // boundary faces from tetgen - - 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(Vector3d(m_V_emb.row(i))); - } - - // add tags - if (!Fs.empty()) { - tag_from_winding_number(); - } else { - assert(!m_img_datas.empty()); - 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 old2new; diff --git a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp index 483a590baf..c437bc3d4f 100644 --- a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp @@ -29,24 +29,6 @@ void delaunay_box_mesh( Vector3d& box_min, Vector3d& box_max); -/** - * @brief Embed a surface in the given volumetric mesh. - * - * @param V_surface Surface vertices. - * @param F_surface Surface faces (triangles). - * @param V_vol Tet vertices. - * @param T_vol Tets. - * @param[out] V_emb Vertices after embedding. - * @param[out] T_emb Tets after embedding. - */ -void embed_surface( - const MatrixXd& V_surface, - const MatrixXi& F_surface, - const MatrixXd& V_vol, - const MatrixXi& T_vol, - MatrixXr& V_emb, - MatrixXi& T_emb); - /** * @brief Embed a surface in the given volumetric mesh. * @@ -67,35 +49,12 @@ void embed_surface( MatrixXi& T_emb, MatrixXi& F_on_surface); -void tag_tets_from_image( - const std::string& filename, - const Matrix4d& xyz2ijk, - const MatrixXd& V, - const MatrixXi& T, - VectorXi& T_tags); - -void tag_tets_from_image( - const ImageData& data, - const Matrix4d& xyz2ijk, - const MatrixXd& V, - const MatrixXi& T, - VectorXi& T_tags); - -void tag_tets_from_images( - const std::vector& data, - const Matrix4d& xyz2ijk, - const MatrixXd& V, - const MatrixXi& T, - MatrixXi& T_tags); - /** * A class for reading an image and converting it into a tet mesh. */ class EmbedSurface { public: - EmbedSurface(const std::vector& img_filenames, const Matrix4d& ijk2xyz); - /** * @brief Input from meshes. */ @@ -119,8 +78,6 @@ class EmbedSurface bool embed_surface(const bool flood_fill = false); - bool embed_surface_tetgen(); - /** * @brief Remove unreferenced vertices. */ @@ -167,8 +124,6 @@ class EmbedSurface private: std::vector m_img_filenames; std::vector m_img_datas; - Matrix4d m_ijk2xyz; // transformation matrix from image to xyz coordinates - Matrix4d m_xyz2ijk; // the surface separating all tags MatrixXd m_V_surface; diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.cpp index 9c7d862934..ad81fe6dd0 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.cpp @@ -631,54 +631,6 @@ void ImageSimulationMesh::write_msh(std::string file, const bool write_envelope) return m_vertex_attribute[i].m_posf; }); - const auto& tets = get_tets(); - msh.add_tets(tets.size(), [&](size_t k) { - auto i = tets[k].tid(*this); - auto vs = oriented_tet_vertices(tets[k]); - std::array data; - for (int j = 0; j < 4; j++) { - data[j] = vs[j].vid(*this); - assert(data[j] < vtx.size()); - } - return data; - }); - - msh.add_tet_vertex_attribute<1>("tv index", [&](size_t i) { - return m_vertex_attribute[i].m_sizing_scalar; - }); - msh.add_tet_attribute<1>("t energy", [&](size_t i) { - return std::cbrt(m_tet_attribute[i].m_quality); - }); - - for (size_t j = 0; j < m_tags_count; ++j) { - msh.add_tet_attribute<1>(fmt::format("tag_{}", j), [&](size_t i) { - return m_tet_attribute[i].tags.count(j) ? 1 : 0; - }); - } - - msh.add_physical_group("ImageVolume"); - - if (m_envelope && write_envelope) { - msh.add_face_vertices(m_V_envelope.size(), [this](size_t k) { return m_V_envelope[k]; }); - msh.add_faces(m_F_envelope.size(), [this](size_t k) { return m_F_envelope[k]; }); - msh.add_physical_group("EnvelopeSurface"); - } - - msh.save(file, true); -} - -void ImageSimulationMesh::write_msh_groups(std::string file, const bool write_envelope) -{ - consolidate_mesh(); - - wmtk::MshData msh; - - const auto& vtx = get_vertices(); - msh.add_tet_vertices(vtx.size(), [&](size_t k) { - auto i = vtx[k].vid(*this); - return m_vertex_attribute[i].m_posf; - }); - const auto& tets = get_tets(); int64_t max_tag = -1; diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h index ffd305d506..979dae0bae 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h @@ -337,11 +337,8 @@ class ImageSimulationMesh : public wmtk::TetMesh ////// Attributes related void write_msh(std::string file, const bool write_envelope = true); - void write_msh_groups(std::string file, const bool write_envelope = true); void output_faces(std::string file, std::function cond); - void init_from_delaunay_box_mesh(const std::vector& vertices); - public: void split_all_edges(); bool split_edge_before(const Tuple& t) override; @@ -484,32 +481,6 @@ class ImageSimulationMesh : public wmtk::TetMesh tbb::enumerable_thread_specific swap_cache; public: - void insertion_by_volumeremesher( - const std::vector& vertices, - const std::vector>& faces, - std::vector& v_rational, - std::vector>& facets_after, - std::vector& is_v_on_input, - std::vector>& tets_after, - std::vector& tet_face_on_input_surface); - - /** - * @brief Init mesh from arrays generated by the volume remesher. - * - * @param v_rational Vertex positions in rational number type. - * @param is_v_on_input A vector of bools identifying if a vertex is on the input. - * @param tets Vector of tet vertex IDs, e.g., [[v0,v1,v2,v3],[...],...]. - * @param tet_face_on_input_surface A vector of bools identifying if the face is on the input. - * The vector contains 4 faces for each tet. If a face is on the input, all incident tets must - * mark that face as input. - * - */ - void init_from_Volumeremesher( - const std::vector& v_rational, - const std::vector& is_v_on_input, - const std::vector>& tets, - const std::vector& tet_face_on_input_surface); - /** * @brief Init from meshes image. * diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index f4d707ea8e..9eef337b0b 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -488,56 +488,6 @@ void ImageSimulationMeshTri::write_msh(std::string file, const bool write_envelo return Vector3d(p2[0], p2[1], 0); }); - const auto faces = get_faces(); - msh.add_faces(faces.size(), [&](size_t k) { - const size_t fid = faces[k].fid(*this); - auto vs = oriented_tri_vertices(faces[k]); - std::array data; - for (size_t j = 0; j < 3; j++) { - data[j] = vs[j].vid(*this); - assert(data[j] < vtx.size()); - } - return data; - }); - - msh.add_face_vertex_attribute<1>("sizing_scalar", [&](size_t i) { - return m_vertex_attribute[i].m_sizing_scalar; - }); - msh.add_face_attribute<1>("quality", [&](size_t i) { return m_face_attribute[i].m_quality; }); - - for (size_t j = 0; j < m_tags_count; ++j) { - msh.add_face_attribute<1>(fmt::format("tag_{}", j), [&](size_t i) { - return m_face_attribute[i].tags.count(j) ? 1 : 0; - }); - } - - msh.add_physical_group("ImageVolume"); - - if (m_envelope && write_envelope) { - msh.add_edge_vertices(m_V_envelope.size(), [this](size_t k) { - return Vector3d(m_V_envelope[k][0], m_V_envelope[k][1], 0); - }); - msh.add_edges(m_E_envelope.size(), [this](size_t k) { return m_E_envelope[k]; }); - msh.add_physical_group("EnvelopeSurface"); - } - - logger().info("Write {}", file); - msh.save(file, true); -} - -void ImageSimulationMeshTri::write_msh_groups(std::string file, const bool write_envelope) -{ - consolidate_mesh(); - - wmtk::MshData msh; - - const auto& vtx = get_vertices(); - msh.add_face_vertices(vtx.size(), [&](size_t k) { - auto i = vtx[k].vid(*this); - Vector2d p2 = m_vertex_attribute[i].m_pos; - return Vector3d(p2[0], p2[1], 0); - }); - const auto& faces = get_faces(); int64_t max_tag = -1; @@ -1491,34 +1441,6 @@ void ImageSimulationMeshTri::smooth_all_vertices(const size_t n_iters) write_vtu(fmt::format("debug_{}", m_debug_print_counter++)); } } - - // re-build envelope - if (m_params.smooth_without_envelope) { - logger().warn("Update envelope"); - MatrixXd V; - V.resize(vert_capacity(), 2); - V.setZero(); - for (size_t i = 0; i < vert_capacity(); ++i) { - const Tuple v = tuple_from_vertex(i); - if (!v.is_valid(*this)) { - continue; - } - const size_t vid = v.vid(*this); - V.row(i) = m_vertex_attribute.at(vid).m_pos; - } - - const auto surf_edges = get_edges_by_condition([](auto& f) { return f.m_is_surface_fs; }); - MatrixXi E; - E.resize(surf_edges.size(), 2); - for (size_t i = 0; i < surf_edges.size(); ++i) { - E.row(i) = Vector2i((int)surf_edges[i][0], (int)surf_edges[i][1]); - } - - m_envelope = nullptr; - m_V_envelope.clear(); - m_E_envelope.clear(); - init_envelope(V, E); - } } bool ImageSimulationMeshTri::smooth_before(const Tuple& t) @@ -1593,7 +1515,7 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) logger().trace("old pos {} -> new pos {}", old_pos.transpose(), VA[vid].m_pos.transpose()); // check surface containment - if (VA[vid].m_is_on_surface && !m_params.smooth_without_envelope) { + if (VA[vid].m_is_on_surface) { // write_vtu_with_energies(fmt::format("debug_smooth_{}", m_debug_print_counter++)); for (size_t i = 1; i < surf_assembles.size(); ++i) { @@ -1656,10 +1578,7 @@ std::shared_ptr ImageSimulationMeshTri::get_envel { const double w = m_s_envelope * m_params.w_envelope; - auto envelope_energy = std::make_shared( - m_envelope_orig, - w, - !m_params.smooth_without_envelope); + auto envelope_energy = std::make_shared(m_envelope_orig, w); return envelope_energy; } diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index eb3aed6e40..39f4e128c2 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp @@ -183,7 +183,6 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool adjust_sizing_field_serial(double max_energy); void write_msh(std::string file, const bool write_envelope = true); - void write_msh_groups(std::string file, const bool write_envelope = true); void write_vtu(const std::string& path) const; void write_vtu_with_energies(const std::string& path) const; diff --git a/components/image_simulation/wmtk/components/image_simulation/Parameters.h b/components/image_simulation/wmtk/components/image_simulation/Parameters.h index 77afedf29b..62f51cf0d7 100644 --- a/components/image_simulation/wmtk/components/image_simulation/Parameters.h +++ b/components/image_simulation/wmtk/components/image_simulation/Parameters.h @@ -1,4 +1,5 @@ #pragma once +#include #include namespace wmtk::components::image_simulation { @@ -18,12 +19,8 @@ struct Parameters // parameters set in `init` function based on mesh bbox double diag_l = -1.; - double diag_l2 = -1.; - double diag_l3 = -1.; - double diag_l4 = -1.; VectorXd box_min; VectorXd box_max; - double vol = -1; // bbox volume double splitting_l2 = -1.; // the lower bound length (squared) for edge split double collapsing_l2 = std::numeric_limits::max(); // the upper bound length (squared) for edge collapse @@ -34,23 +31,55 @@ struct Parameters bool debug_output = false; bool perform_sanity_checks = false; - bool smooth_without_envelope = false; - // weighting terms for the optimization double w_amips = 1e-4; double w_envelope = 0; std::string operation = "remeshing"; + bool skip_simplify = true; + bool use_sample_envelope = false; + int NUM_THREADS = 0; + int max_its = 80; + bool write_vtu = false; + bool write_envelope = true; + + Parameters() = default; + + Parameters(const nlohmann::json& json_params) + { + output_path = json_params["output"]; + skip_simplify = json_params["skip_simplify"]; + use_sample_envelope = json_params["use_sample_envelope"]; + NUM_THREADS = json_params["num_threads"]; + max_its = json_params["max_iterations"]; + write_vtu = json_params["write_vtu"]; + write_envelope = json_params["write_envelope"]; + + epsr = json_params["eps_rel"]; + eps = json_params["eps"]; + lr = json_params["length_rel"]; + l = json_params["length"]; + stop_energy = json_params["stop_energy"]; + stop_at_float = json_params["stop_at_float"]; + preserve_topology = json_params["preserve_topology"]; + + epsr_simplify = json_params["eps_simplify_rel"]; + eps_simplify = json_params["eps_simplify"]; + + w_amips = json_params["w_amips"]; + + debug_output = json_params["DEBUG_output"]; + perform_sanity_checks = json_params["DEBUG_sanity_checks"]; + + operation = json_params["operation"]; + } + void init(const VectorXd& min_, const VectorXd& max_) { box_min = min_; box_max = max_; - vol = (box_max - box_min).prod(); diag_l = (box_max - box_min).norm(); - diag_l2 = diag_l * diag_l; - diag_l3 = diag_l * diag_l * diag_l; - diag_l4 = diag_l * diag_l * diag_l * diag_l; if (l > 0) lr = l / diag_l; else diff --git a/components/image_simulation/wmtk/components/image_simulation/Smooth.cpp b/components/image_simulation/wmtk/components/image_simulation/Smooth.cpp index 7038f7e263..3a9f07946a 100644 --- a/components/image_simulation/wmtk/components/image_simulation/Smooth.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/Smooth.cpp @@ -111,7 +111,7 @@ bool ImageSimulationMesh::smooth_after(const Tuple& t) logger().trace("old pos {} -> new pos {}", old_pos, VA[vid].m_posf); // check surface containment - if (VA[vid].m_is_on_surface && !m_params.smooth_without_envelope) { + if (VA[vid].m_is_on_surface) { // write_vtu_with_energies(fmt::format("debug_smooth_{}", debug_print_counter++)); const simplex::SimplexCollection surf_assembles = get_surface_faces_for_vertex(vid); for (size_t i = 0; i < surf_assembles.faces().size(); ++i) { @@ -185,35 +185,6 @@ void ImageSimulationMesh::smooth_all_vertices(const size_t n_iters = 1) write_vtu(fmt::format("debug_{}", m_debug_print_counter++)); } } - - // re-build envelope - if (m_params.smooth_without_envelope) { - logger().warn("Update envelope"); - MatrixXd V; - V.resize(vert_capacity(), 3); - V.setZero(); - for (size_t i = 0; i < vert_capacity(); ++i) { - const Tuple v = tuple_from_vertex(i); - if (!v.is_valid(*this)) { - continue; - } - const size_t vid = v.vid(*this); - V.row(i) = m_vertex_attribute.at(vid).m_posf; - } - - const auto surf_faces = get_faces_by_condition([](auto& f) { return f.m_is_surface_fs; }); - MatrixXi F; - F.resize(surf_faces.size(), 3); - for (size_t i = 0; i < surf_faces.size(); ++i) { - F.row(i) = - Vector3i((int)surf_faces[i][0], (int)surf_faces[i][1], (int)surf_faces[i][2]); - } - - m_envelope = nullptr; - m_V_envelope.clear(); - m_F_envelope.clear(); - init_envelope(V, F); - } } std::vector> ImageSimulationMesh::get_amips_assembles(const Tuple& t) const @@ -271,8 +242,7 @@ std::shared_ptr ImageSimulationMesh::get_envelope m_vertex_attribute[vid].m_order); } - auto envelope_energy = - std::make_shared(env, w, !m_params.smooth_without_envelope); + auto envelope_energy = std::make_shared(env, w); return envelope_energy; } diff --git a/components/image_simulation/wmtk/components/image_simulation/VolumemesherInsertion.cpp b/components/image_simulation/wmtk/components/image_simulation/VolumemesherInsertion.cpp index 11bdaae4ae..0ee5321c38 100644 --- a/components/image_simulation/wmtk/components/image_simulation/VolumemesherInsertion.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/VolumemesherInsertion.cpp @@ -8,708 +8,6 @@ namespace wmtk::components::image_simulation { -std::vector> ImageSimulationMesh::triangulate_polygon_face( - std::vector points) -{ - // triangulate weak convex polygons - std::vector> triangulated_faces; - - std::vector> points_vector; - for (int i = 0; i < points.size(); i++) { - points_vector.push_back(std::pair(points[i], i)); - } - - // find the first colinear ABC with nonlinear BCD and delete C from vector - while (points_vector.size() > 3) { - bool no_colinear = true; - for (int i = 0; i < points_vector.size(); i++) { - auto cur = points_vector[i]; - auto next = points_vector[(i + 1) % points_vector.size()]; - auto prev = points_vector[(i + points_vector.size() - 1) % points_vector.size()]; - auto nextnext = points_vector[(i + 2) % points_vector.size()]; - - Vector3r a = cur.first - prev.first; - Vector3r b = next.first - cur.first; - Vector3r c = nextnext.first - next.first; - - if (((a[0] * b[1] - a[1] * b[0]) == 0 && (a[1] * b[2] - a[2] * b[1]) == 0 && - (a[0] * b[2] - a[2] * b[0]) == 0) && - (((b[0] * c[1] - b[1] * c[0]) != 0 || (b[1] * c[2] - b[2] * c[1]) != 0 || - (b[0] * c[2] - b[2] * c[0]) != 0))) { - no_colinear = false; - std::array t = { - size_t(cur.second), - size_t(next.second), - size_t(nextnext.second)}; - triangulated_faces.push_back(t); - points_vector.erase(points_vector.begin() + ((i + 1) % points_vector.size())); - break; - } else { - continue; - } - } - - if (no_colinear) break; - } - - // cleanup convex polygon - while (points_vector.size() >= 3) { - std::array t = { - size_t(points_vector[0].second), - size_t(points_vector[1].second), - size_t(points_vector[points_vector.size() - 1].second)}; - triangulated_faces.push_back(t); - points_vector.erase(points_vector.begin()); - } - - return triangulated_faces; -} - -// we have the vertices and triangles -// to generate what we need for volumemesher -// coords, ncoords, tri_idx, ntri_idx - -// embed input surface on generated back ground mesh - -void ImageSimulationMesh::init_from_delaunay_box_mesh(const std::vector& vertices) -{ - ///points for delaunay - std::vector points(vertices.size()); - // add points from surface - for (int i = 0; i < vertices.size(); i++) { - for (int j = 0; j < 3; j++) points[i][j] = vertices[i][j]; - } - - // bbox - double delta = m_params.diag_l / 15.0; - Vector3d box_min( - m_params.box_min[0] - delta, - m_params.box_min[1] - delta, - m_params.box_min[2] - delta); - Vector3d box_max( - m_params.box_max[0] + delta, - m_params.box_max[1] + delta, - m_params.box_max[2] + delta); - - // add corners of domain - for (int i = 0; i < 8; i++) { - Vector3d p; - std::bitset a(i); - for (int j = 0; j < 3; j++) { - if (a.test(j)) { - p[j] = box_max[j]; - } else { - p[j] = box_min[j]; - } - } - points.push_back({{p[0], p[1], p[2]}}); - } - - const double voxel_resolution = m_params.diag_l / 20.0; - std::array N; // number of grid points per dimension - std::array 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, 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 = m_envelope->nearest_point(p, n); - - if (sqd < min_dis) { - continue; - } - points.push_back({{ds[0][i], ds[1][j], ds[2][k]}}); - } - } - } - - m_params.box_min = box_min; - m_params.box_max = box_max; - - ///delaunay - auto [unused_points, tets] = delaunay::delaunay3D(points); - logger().info("after delauney tets.size() {} points.size() {}", tets.size(), points.size()); - - // conn - init(points.size(), tets); - // attr - m_vertex_attribute.m_attributes.resize(points.size()); - m_tet_attribute.m_attributes.resize(tets.size()); - m_face_attribute.m_attributes.resize(tets.size() * 4); - for (int i = 0; i < vert_capacity(); i++) { - m_vertex_attribute[i].m_pos = Vector3r(points[i][0], points[i][1], points[i][2]); - m_vertex_attribute[i].m_posf = Vector3d(points[i][0], points[i][1], points[i][2]); - } -} - -void ImageSimulationMesh::insertion_by_volumeremesher( - const std::vector& vertices, - const std::vector>& faces, - std::vector& v_rational, - std::vector>& facets_after, - std::vector& is_v_on_input, - std::vector>& tets_after, - std::vector& tet_face_on_input_surface) -{ - std::cout << "vertices size: " << vertices.size() << std::endl; - std::cout << "faces size: " << faces.size() << std::endl; - - // generate background mesh - init_from_delaunay_box_mesh(vertices); - - // prepare tet vertices and tet index info - - auto tet_vers = get_vertices(); - auto tets = get_tets(); - std::vector tet_ver_coord(3 * tet_vers.size()); - std::vector tet_index(4 * tets.size()); - std::cout << "tetver size: " << tet_vers.size() << std::endl; - std::cout << "tet size: " << tets.size() << std::endl; - - for (int i = 0; i < tet_vers.size(); ++i) { - tet_ver_coord[3 * i] = m_vertex_attribute[i].m_posf[0]; - tet_ver_coord[3 * i + 1] = m_vertex_attribute[i].m_posf[1]; - tet_ver_coord[3 * i + 2] = m_vertex_attribute[i].m_posf[2]; - } - - for (int i = 0; i < tets.size(); ++i) { - auto tet_vids = oriented_tet_vids(tets[i]); - tet_index[4 * i] = tet_vids[0]; - tet_index[4 * i + 1] = tet_vids[1]; - tet_index[4 * i + 2] = tet_vids[2]; - tet_index[4 * i + 3] = tet_vids[3]; - } - - // prepare input surfaces info - std::vector tri_ver_coord(3 * vertices.size()); - std::vector tri_index(3 * faces.size()); - - for (int i = 0; i < vertices.size(); ++i) { - tri_ver_coord[3 * i] = vertices[i][0]; - tri_ver_coord[3 * i + 1] = vertices[i][1]; - tri_ver_coord[3 * i + 2] = vertices[i][2]; - } - - for (int i = 0; i < faces.size(); ++i) { - tri_index[3 * i] = faces[i][0]; - tri_index[3 * i + 1] = faces[i][1]; - tri_index[3 * i + 2] = faces[i][2]; - } - - std::cout << tri_ver_coord.size() << std::endl; - std::cout << tri_index.size() << std::endl; - std::cout << tet_ver_coord.size() << std::endl; - std::cout << tet_index.size() << std::endl; - - std::vector embedded_vertices; - std::vector embedded_facets; - std::vector embedded_cells; - std::vector embedded_facets_on_input; - - std::vector> out_tets; - std::vector final_tets_parent; - std::vector cells_with_faces_on_input; - std::vector> final_tets_parent_faces; - - { - logger().info("Check degenerate before embedding"); - for (int i = 0; i < tri_index.size(); i += 3) { - int id0 = tri_index[i + 0]; - int id1 = tri_index[i + 1]; - int id2 = tri_index[i + 2]; - Eigen::Vector3d v0; - Eigen::Vector3d v1; - Eigen::Vector3d v2; - v0 << tri_ver_coord[3 * id0 + 0], tri_ver_coord[3 * id0 + 1], - tri_ver_coord[3 * id0 + 2]; - v1 << tri_ver_coord[3 * id1 + 0], tri_ver_coord[3 * id1 + 1], - tri_ver_coord[3 * id1 + 2]; - v2 << tri_ver_coord[3 * id2 + 0], tri_ver_coord[3 * id2 + 1], - tri_ver_coord[3 * id2 + 2]; - - if (utils::predicates::is_degenerate(v0, v1, v2)) { - logger().warn( - "Face ({}, {}, {}) is collinear!", - v0.transpose(), - v1.transpose(), - v2.transpose()); - } - } - - for (int i = 0; i < tet_index.size(); i += 4) { - int id0 = tet_index[i + 0]; - int id1 = tet_index[i + 1]; - int id2 = tet_index[i + 2]; - int id3 = tet_index[i + 3]; - Eigen::Vector3d v0; - Eigen::Vector3d v1; - Eigen::Vector3d v2; - Eigen::Vector3d v3; - v0 << tet_ver_coord[3 * id0 + 0], tet_ver_coord[3 * id0 + 1], - tet_ver_coord[3 * id0 + 2]; - v1 << tet_ver_coord[3 * id1 + 0], tet_ver_coord[3 * id1 + 1], - tet_ver_coord[3 * id1 + 2]; - v2 << tet_ver_coord[3 * id2 + 0], tet_ver_coord[3 * id2 + 1], - tet_ver_coord[3 * id2 + 2]; - v3 << tet_ver_coord[3 * id3 + 0], tet_ver_coord[3 * id3 + 1], - tet_ver_coord[3 * id3 + 2]; - - if (utils::predicates::is_degenerate(v0, v1, v2, v3)) { - logger().warn( - "Tet ({}, {}, {}) is coplanar!", - v0.transpose(), - v1.transpose(), - v2.transpose(), - v3.transpose()); - } - } - logger().info("done"); - } - - // volumeremesher embed - vol_rem::embed_tri_in_poly_mesh( - tri_ver_coord, - tri_index, - tet_ver_coord, - tet_index, - embedded_vertices, - embedded_facets, - embedded_cells, - out_tets, - final_tets_parent, - embedded_facets_on_input, - cells_with_faces_on_input, - final_tets_parent_faces, - true); - - //v_rational.resize(embedded_vertices.size() / 3); - for (int i = 0; i < embedded_vertices.size() / 3; i++) { - v_rational.push_back(Vector3r()); -#ifdef USE_GNU_GMP_CLASSES - v_rational.back()[0].init(embedded_vertices[3 * i].get_mpq_t()); - v_rational.back()[1].init(embedded_vertices[3 * i + 1].get_mpq_t()); - v_rational.back()[2].init(embedded_vertices[3 * i + 2].get_mpq_t()); -#else - v_rational.back()[0].init_from_bin(embedded_vertices[3 * i].get_str()); - v_rational.back()[1].init_from_bin(embedded_vertices[3 * i + 1].get_str()); - v_rational.back()[2].init_from_bin(embedded_vertices[3 * i + 2].get_str()); -#endif - } - - std::vector> polygon_faces; - int polycnt = 0; - for (int i = 0; i < embedded_facets.size(); i++) { - int polysize = embedded_facets[i]; - std::vector polygon; - for (int j = i + 1; j <= i + polysize; j++) { - polygon.push_back(embedded_facets[j]); - } - polycnt++; - polygon_faces.push_back(polygon); - i += polysize; - } - - // test polygon faces - only the triangles - for (int i = 0; i < polygon_faces.size(); ++i) { - if (polygon_faces[i].size() != 3) { - continue; - } - - const Vector3r v0 = v_rational[polygon_faces[i][0]]; - const Vector3r v1 = v_rational[polygon_faces[i][1]]; - const Vector3r v2 = v_rational[polygon_faces[i][2]]; - const Vector3r v01 = (v1 - v0); - const Vector3r v02 = (v2 - v0); - const Vector3r r = v01.cross(v02); - if (r[0] == 0 && r[1] == 0 && r[2] == 0) { - logger().error("Collinear triangle in polygon faces. ID = {}", i); - logger().error("v0 = {}", to_double(v0)); - logger().error("v1 = {}", to_double(v1)); - logger().error("v2 = {}", to_double(v2)); - } - } - - std::vector> polygon_cells; - std::vector> tets_final; - for (int i = 0; i < embedded_cells.size(); i++) { - std::vector polygon_cell; - int cellsize = embedded_cells[i]; - for (int j = i + 1; j <= i + cellsize; j++) { - polygon_cell.push_back(embedded_cells[j]); - } - polygon_cells.push_back(polygon_cell); - i += cellsize; - } - - std::cout << "polygon cells num: " << polygon_cells.size() << std::endl; - - std::vector polygon_faces_on_input_surface(polygon_faces.size(), false); - - for (int i = 0; i < embedded_facets_on_input.size(); i++) { - polygon_faces_on_input_surface[embedded_facets_on_input[i]] = true; - } - - std::vector> triangulated_faces; - std::vector triangulated_faces_on_input; - std::vector> map_poly_to_tri_face(polygon_faces.size()); - - int poly_cnt = 0; - - // triangulate polygon faces - for (int i = 0; i < polygon_faces.size(); i++) { - // already clipped in other polygon - if (map_poly_to_tri_face[i].size() != 0) continue; - - // new polygon face to clip - std::vector> clipped_indices; - std::vector poly_coordinates; - std::vector polygon_face = polygon_faces[i]; - assert(polygon_face.size() >= 3); - - if (polygon_face.size() == 3) { - // already a triangle - std::array triangle_face = { - polygon_face[0], - polygon_face[1], - polygon_face[2]}; - int idx = triangulated_faces.size(); - triangulated_faces.push_back(triangle_face); - if (polygon_faces_on_input_surface[i]) { - triangulated_faces_on_input.push_back(true); - } else { - triangulated_faces_on_input.push_back(false); - } - map_poly_to_tri_face[i].push_back(idx); - } else { - poly_cnt++; - for (int j = 0; j < polygon_faces[i].size(); j++) { - poly_coordinates.push_back(v_rational[polygon_face[j]]); - } - - clipped_indices = triangulate_polygon_face(poly_coordinates); - for (int j = 0; j < clipped_indices.size(); j++) { - // need to map oldface index to new face indices - std::array triangle_face = { - polygon_face[clipped_indices[j][0]], - polygon_face[clipped_indices[j][1]], - polygon_face[clipped_indices[j][2]]}; - int idx = triangulated_faces.size(); - triangulated_faces.push_back(triangle_face); - - // track input faces - if (polygon_faces_on_input_surface[i]) { - triangulated_faces_on_input.push_back(true); - } else { - triangulated_faces_on_input.push_back(false); - } - map_poly_to_tri_face[i].push_back(idx); - } - } - } - - std::cout << "poly_cnt:" << poly_cnt << std::endl; - std::cout << "finish triangulation" << std::endl; - std::cout << "vertice before tetra num: " << v_rational.size() << std::endl; - - int was_tet_cnt = 0; - for (int i = 0; i < polygon_cells.size(); i++) { - auto polygon_cell = polygon_cells[i]; - - // get polygon vertices - std::vector polygon_vertices; - for (auto f : polygon_cell) { - for (auto v : polygon_faces[f]) { - polygon_vertices.push_back(v); - } - } - wmtk::vector_unique(polygon_vertices); - - // compute number of triangle faces - int num_faces = 0; - for (auto f : polygon_cell) { - num_faces += map_poly_to_tri_face[f].size(); - } - - // polygon already a tet - if (num_faces == 4) { - was_tet_cnt++; - assert(polygon_vertices.size() == 4); - // get the correct orientation here - size_t v0 = polygon_faces[polygon_cell[0]][0]; - size_t v1 = polygon_faces[polygon_cell[0]][1]; - size_t v2 = polygon_faces[polygon_cell[0]][2]; - size_t v3; - for (auto v : polygon_faces[polygon_cell[1]]) { - if (v != v0 && v != v1 && v != v2) { - v3 = v; - break; - } - } - - std::array tetra = {v0, v1, v2, v3}; - - // if inverted then fix the orientation - Vector3r v0v1 = v_rational[v1] - v_rational[v0]; - Vector3r v0v2 = v_rational[v2] - v_rational[v0]; - Vector3r v0v3 = v_rational[v3] - v_rational[v0]; - if ((v0v1.cross(v0v2)).dot(v0v3) < 0) { - tetra = {v1, v0, v2, v3}; - } - - // push the tet to final queue; - tets_final.push_back(tetra); - - std::set local_f1 = {tetra[0], tetra[1], tetra[2]}; - std::set local_f2 = {tetra[0], tetra[2], tetra[3]}; - std::set local_f3 = {tetra[0], tetra[1], tetra[3]}; - std::set local_f4 = {tetra[1], tetra[2], tetra[3]}; - - // track surface need to be fixed - bool tet_face_on_input[4]; - for (auto f : polygon_cell) { - std::set f_vs = { - polygon_faces[f][0], - polygon_faces[f][1], - polygon_faces[f][2]}; - - int local_f_idx; - - // decide which face it is - - if (f_vs == local_f1) { - local_f_idx = 0; - } else if (f_vs == local_f2) { - local_f_idx = 1; - } else if (f_vs == local_f3) { - local_f_idx = 2; - } else { - local_f_idx = 3; - } - - tet_face_on_input[local_f_idx] = polygon_faces_on_input_surface[f]; - } - - for (int k = 0; k < 4; k++) { - tet_face_on_input_surface.push_back(tet_face_on_input[k]); - } - continue; - } - - // compute centroid - Vector3r centroid(0, 0, 0); - for (auto v : polygon_vertices) { - centroid = centroid + v_rational[v]; - } - centroid = centroid / polygon_vertices.size(); - - // trahedralize - size_t centroid_idx = v_rational.size(); - v_rational.push_back(centroid); - - for (auto f : polygon_cell) { - for (auto t : map_poly_to_tri_face[f]) { - std::array tetra = { - triangulated_faces[t][0], - triangulated_faces[t][1], - triangulated_faces[t][2], - centroid_idx}; - // check inverted tet and fix - Vector3r v0v1 = v_rational[tetra[1]] - v_rational[tetra[0]]; - Vector3r v0v2 = v_rational[tetra[2]] - v_rational[tetra[0]]; - Vector3r v0v3 = v_rational[tetra[3]] - v_rational[tetra[0]]; - if ((v0v1.cross(v0v2)).dot(v0v3) < 0) { - tetra = { - triangulated_faces[t][1], - triangulated_faces[t][0], - triangulated_faces[t][2], - centroid_idx}; - } - - tets_final.push_back(tetra); - tet_face_on_input_surface.push_back(triangulated_faces_on_input[t]); - tet_face_on_input_surface.push_back(false); - tet_face_on_input_surface.push_back(false); - tet_face_on_input_surface.push_back(false); - } - } - } - - std::cout << "polygon was tet num: " << was_tet_cnt << std::endl; - std::cout << "vertices final num: " << v_rational.size() << std::endl; - std::cout << "tets final num: " << tets_final.size() << std::endl; - - std::cout << "track face size: " << tet_face_on_input_surface.size() << std::endl; - - facets_after = triangulated_faces; - tets_after = tets_final; - - // track vertices on input - is_v_on_input.reserve(v_rational.size()); - for (int i = 0; i < v_rational.size(); i++) { - is_v_on_input.push_back(false); - } - for (int i = 0; i < triangulated_faces.size(); i++) { - if (triangulated_faces_on_input[i]) { - is_v_on_input[triangulated_faces[i][0]] = true; - is_v_on_input[triangulated_faces[i][1]] = true; - is_v_on_input[triangulated_faces[i][2]] = true; - } - } - - size_t on_surface_v_cnt = 0; - for (size_t i = 0; i < is_v_on_input.size(); i++) { - if (is_v_on_input[i]) on_surface_v_cnt++; - } - - for (const auto& t : tets_after) { - const Vector3r p0 = v_rational[t[0]]; - const Vector3r p1 = v_rational[t[1]]; - const Vector3r p2 = v_rational[t[2]]; - const Vector3r p3 = v_rational[t[3]]; - Vector3r n = (p1 - p0).cross(p2 - p0); - Vector3r d = p3 - p0; - auto res = n.dot(d); - if (res <= 0) { - logger().error("Inverted tet {}", t); - } - } - std::cout << "v on surface vector size: " << is_v_on_input.size(); - std::cout << "v on surface: " << on_surface_v_cnt << std::endl; -} - -void ImageSimulationMesh::init_from_Volumeremesher( - const std::vector& v_rational, - const std::vector& is_v_on_input, - const std::vector>& tets, - const std::vector& tet_face_on_input_surface) -{ - assert(tet_face_on_input_surface.size() == 4 * tets.size()); - - init_with_isolated_vertices(v_rational.size(), tets); - assert(check_mesh_connectivity_validity()); - - m_vertex_attribute.m_attributes.resize(v_rational.size()); - m_tet_attribute.m_attributes.resize(tets.size()); - m_face_attribute.m_attributes.resize(tets.size() * 4); - - for (int i = 0; i < vert_capacity(); i++) { - m_vertex_attribute[i].m_pos = v_rational[i]; - m_vertex_attribute[i].m_posf = to_double(v_rational[i]); - } - - // check here - for (size_t i = 0; i < tet_face_on_input_surface.size(); i++) { - if (tet_face_on_input_surface[i]) { - m_face_attribute[i].m_is_surface_fs = 1; - } - } - - const auto faces = get_faces(); - std::cout << "faces size: " << faces.size() << std::endl; - for (const Tuple& f : faces) { - SmartTuple ff(*this, f); - const size_t fid = ff.fid(); - const size_t v1 = ff.vid(); - const size_t v2 = ff.switch_vertex().vid(); - const size_t v3 = ff.switch_edge().switch_vertex().vid(); - if (m_face_attribute[fid].m_is_surface_fs == 1) { - assert(is_v_on_input[v1] && is_v_on_input[v2] && is_v_on_input[v3]); - m_vertex_attribute[v1].m_is_on_surface = true; - m_vertex_attribute[v2].m_is_on_surface = true; - m_vertex_attribute[v3].m_is_on_surface = true; - } - } - - // track bounding box - for (size_t i = 0; i < faces.size(); i++) { - const auto vs = get_face_vertices(faces[i]); - std::array vids = {{vs[0].vid(*this), vs[1].vid(*this), vs[2].vid(*this)}}; - int on_bbox = -1; - for (int k = 0; k < 3; k++) { - if (m_vertex_attribute[vids[0]].m_pos[k] == m_params.box_min[k] && - m_vertex_attribute[vids[1]].m_pos[k] == m_params.box_min[k] && - m_vertex_attribute[vids[2]].m_pos[k] == m_params.box_min[k]) { - on_bbox = k * 2; - break; - } - if (m_vertex_attribute[vids[0]].m_pos[k] == m_params.box_max[k] && - m_vertex_attribute[vids[1]].m_pos[k] == m_params.box_max[k] && - m_vertex_attribute[vids[2]].m_pos[k] == m_params.box_max[k]) { - on_bbox = k * 2 + 1; - break; - } - } - if (on_bbox < 0) { - continue; - } - assert(!faces[i].switch_tetrahedron(*this)); // face must be on boundary - - const size_t fid = faces[i].fid(*this); - m_face_attribute[fid].m_is_bbox_fs = on_bbox; - - for (const size_t vid : vids) { - m_vertex_attribute[vid].on_bbox_faces.push_back(on_bbox); - } - } - - for_each_vertex( - [&](auto& v) { wmtk::vector_unique(m_vertex_attribute[v.vid(*this)].on_bbox_faces); }); - - // track open boundaries - find_order_2_edges(); - - int open_boundary_cnt = 0; - for (const Tuple& e : get_edges()) { - if (is_order_2_edge(e)) { - open_boundary_cnt++; - } - } - wmtk::logger().info("#open boundary edges: {}", open_boundary_cnt); - - // // rounding - size_t cnt_round = 0; - - const auto vertices = get_vertices(); - for (const Tuple& v : vertices) { - if (round(v)) { - cnt_round++; - } - } - - if (cnt_round != vertices.size()) { - log_and_throw_error( - "Could not round all vertices in tet mesh. Rounded {}/{}", - cnt_round, - vertices.size()); - } - - // init qualities - for_each_tetra( - [this](const Tuple& t) { m_tet_attribute[t.tid(*this)].m_quality = get_quality(t); }); -} - void ImageSimulationMesh::init_from_image( const MatrixXr& V, const MatrixXi& T, @@ -985,7 +283,6 @@ void ImageSimulationMesh::init_surfaces_and_boundaries() } } - void ImageSimulationMesh::find_order_2_edges() { // for open boundary envelope diff --git a/components/image_simulation/wmtk/components/image_simulation/extract_soup.cpp b/components/image_simulation/wmtk/components/image_simulation/extract_soup.cpp deleted file mode 100644 index 305c207334..0000000000 --- a/components/image_simulation/wmtk/components/image_simulation/extract_soup.cpp +++ /dev/null @@ -1,395 +0,0 @@ -#include "extract_soup.hpp" - -#include -//#include -#include -#include -#include - -#include "ImageSimulationMesh.h" - -namespace { - -void read_array_data( - std::vector>>& data, - const std::string& filename) -{ - std::ifstream file(filename, std::ios::binary); - if (!file.is_open()) { - std::cerr << "Error: Unable to open file " << filename << " for reading." << std::endl; - return; - } - - int dim1, dim2, dim3; - file.read(reinterpret_cast(&dim1), sizeof(int)); - file.read(reinterpret_cast(&dim2), sizeof(int)); - file.read(reinterpret_cast(&dim3), sizeof(int)); - - data.resize(dim3, std::vector>(dim2, std::vector(dim1))); - - for (int k = 0; k < dim1; ++k) { - for (int j = 0; j < dim2; ++j) { - for (int i = 0; i < dim3; ++i) { - size_t value; - file.read(reinterpret_cast(&value), sizeof(size_t)); - data[i][j][k] = value; - } - } - } -} - -} // namespace - - -namespace wmtk::components::image_simulation { - -void read_array_data_ascii( - std::vector>>& data, - const std::string& filename) -{ - std::ifstream file(filename, std::ios::binary); - if (!file.is_open()) { - wmtk::log_and_throw_error("Error: Unable to open file {}", filename); - } - - int dim1, dim2, dim3; - file >> dim1; - file >> dim2; - file >> dim3; - - data.resize(dim3, std::vector>(dim2, std::vector(dim1))); - - for (int k = 0; k < dim1; ++k) { - for (int j = 0; j < dim2; ++j) { - for (int i = 0; i < dim3; ++i) { - size_t value; - file >> value; - data[i][j][k] = value; - } - } - } -} - -void extract_triangle_soup_from_image(std::string filename, Eigen::MatrixXi& F, Eigen::MatrixXd& V) -{ - // std::vector>> data; - // read_array_data(data, filename); - std::vector>> data; - read_array_data_ascii(data, filename); - - MatrixXi F_tags; - extract_triangle_soup_from_image(data, F, V, F_tags); -} - -void extract_triangle_soup_from_image( - const std::vector>>& data, - MatrixXi& F, - MatrixXd& V, - MatrixXi& F_tags) -{ - // get the number of triangles - size_t tri_num = 0; - for (size_t i = 0; i < data.size() - 1; i++) { - for (size_t j = 0; j < data[0].size() - 1; j++) { - for (size_t k = 0; k < data[0][0].size() - 1; k++) { - size_t cur_data = data[i][j][k]; - size_t right_data = data[i][j][k + 1]; - size_t ahead_data = data[i][j + 1][k]; - size_t top_data = data[i + 1][j][k]; - if (cur_data != right_data) { - tri_num += 2; - } - if (cur_data != ahead_data) { - tri_num += 2; - } - if (cur_data != top_data) { - tri_num += 2; - } - } - } - } - - F.resize(tri_num, 3); - V.resize(3 * tri_num, 3); - F_tags.resize(tri_num, 2); - - // build triangles whenever the data changes between two cells - size_t cur_v_id = 0; - size_t cur_f_id = 0; - for (size_t i = 0; i < data.size() - 1; i++) { - for (size_t j = 0; j < data[0].size() - 1; j++) { - for (size_t k = 0; k < data[0][0].size() - 1; k++) { - size_t cur_data = data[i][j][k]; - size_t right_data = data[i][j][k + 1]; - size_t ahead_data = data[i][j + 1][k]; - size_t top_data = data[i + 1][j][k]; - Eigen::RowVector3d v0(i, (j + 1), k); - Eigen::RowVector3d v1(i, (j + 1), (k + 1)); - Eigen::RowVector3d v2(i, j, (k + 1)); - Eigen::RowVector3d v3((i + 1), (j + 1), k); - Eigen::RowVector3d v4((i + 1), (j + 1), (k + 1)); - Eigen::RowVector3d v5((i + 1), j, (k + 1)); - Eigen::RowVector3d v6((i + 1), j, k); - - if (cur_data > right_data) { - size_t vid0, vid1, vid2, vid3, vid4, vid5, fid0, fid1; - vid0 = cur_v_id; - vid1 = vid0 + 1; - vid2 = vid1 + 1; - vid3 = vid2 + 1; - vid4 = vid3 + 1; - vid5 = vid4 + 1; - fid0 = cur_f_id; - fid1 = fid0 + 1; - cur_v_id += 6; - cur_f_id += 2; - F.row(fid0) << vid0, vid1, vid2; - F.row(fid1) << vid3, vid4, vid5; - V.row(vid0) = v2; - V.row(vid1) = v5; - V.row(vid2) = v4; - V.row(vid3) = v2; - V.row(vid4) = v4; - V.row(vid5) = v1; - F_tags.row(fid0) = Vector2i(right_data, cur_data); - F_tags.row(fid1) = Vector2i(right_data, cur_data); - } else if (cur_data < right_data) { - size_t vid0, vid1, vid2, vid3, vid4, vid5, fid0, fid1; - vid0 = cur_v_id; - vid1 = vid0 + 1; - vid2 = vid1 + 1; - vid3 = vid2 + 1; - vid4 = vid3 + 1; - vid5 = vid4 + 1; - fid0 = cur_f_id; - fid1 = fid0 + 1; - cur_v_id += 6; - cur_f_id += 2; - F.row(fid0) << vid0, vid1, vid2; - F.row(fid1) << vid3, vid4, vid5; - V.row(vid0) = v1; - V.row(vid1) = v4; - V.row(vid2) = v2; - V.row(vid3) = v4; - V.row(vid4) = v5; - V.row(vid5) = v2; - F_tags.row(fid0) = Vector2i(cur_data, right_data); - F_tags.row(fid1) = Vector2i(cur_data, right_data); - } - - if (cur_data > ahead_data) { - size_t vid0, vid1, vid2, vid3, vid4, vid5, fid0, fid1; - vid0 = cur_v_id; - vid1 = vid0 + 1; - vid2 = vid1 + 1; - vid3 = vid2 + 1; - vid4 = vid3 + 1; - vid5 = vid4 + 1; - fid0 = cur_f_id; - fid1 = fid0 + 1; - cur_v_id += 6; - cur_f_id += 2; - F.row(fid0) << vid0, vid1, vid2; - F.row(fid1) << vid3, vid4, vid5; - V.row(vid0) = v0; - V.row(vid1) = v1; - V.row(vid2) = v3; - V.row(vid3) = v1; - V.row(vid4) = v4; - V.row(vid5) = v3; - F_tags.row(fid0) = Vector2i(ahead_data, cur_data); - F_tags.row(fid1) = Vector2i(ahead_data, cur_data); - } else if (cur_data < ahead_data) { - size_t vid0, vid1, vid2, vid3, vid4, vid5, fid0, fid1; - vid0 = cur_v_id; - vid1 = vid0 + 1; - vid2 = vid1 + 1; - vid3 = vid2 + 1; - vid4 = vid3 + 1; - vid5 = vid4 + 1; - fid0 = cur_f_id; - fid1 = fid0 + 1; - cur_v_id += 6; - cur_f_id += 2; - F.row(fid0) << vid0, vid1, vid2; - F.row(fid1) << vid3, vid4, vid5; - V.row(vid0) = v0; - V.row(vid1) = v3; - V.row(vid2) = v1; - V.row(vid3) = v3; - V.row(vid4) = v4; - V.row(vid5) = v1; - F_tags.row(fid0) = Vector2i(cur_data, ahead_data); - F_tags.row(fid1) = Vector2i(cur_data, ahead_data); - } - - if (cur_data > top_data) { - size_t vid0, vid1, vid2, vid3, vid4, vid5, fid0, fid1; - vid0 = cur_v_id; - vid1 = vid0 + 1; - vid2 = vid1 + 1; - vid3 = vid2 + 1; - vid4 = vid3 + 1; - vid5 = vid4 + 1; - fid0 = cur_f_id; - fid1 = fid0 + 1; - cur_v_id += 6; - cur_f_id += 2; - F.row(fid0) << vid0, vid1, vid2; - F.row(fid1) << vid3, vid4, vid5; - V.row(vid0) = v6; - V.row(vid1) = v3; - V.row(vid2) = v5; - V.row(vid3) = v3; - V.row(vid4) = v4; - V.row(vid5) = v5; - F_tags.row(fid0) = Vector2i(top_data, cur_data); - F_tags.row(fid1) = Vector2i(top_data, cur_data); - } else if (cur_data < top_data) { - size_t vid0, vid1, vid2, vid3, vid4, vid5, fid0, fid1; - vid0 = cur_v_id; - vid1 = vid0 + 1; - vid2 = vid1 + 1; - vid3 = vid2 + 1; - vid4 = vid3 + 1; - vid5 = vid4 + 1; - fid0 = cur_f_id; - fid1 = fid0 + 1; - cur_v_id += 6; - cur_f_id += 2; - F.row(fid0) << vid0, vid1, vid2; - F.row(fid1) << vid3, vid4, vid5; - V.row(vid0) = v3; - V.row(vid1) = v5; - V.row(vid2) = v4; - V.row(vid3) = v3; - V.row(vid4) = v6; - V.row(vid5) = v5; - F_tags.row(fid0) = Vector2i(cur_data, top_data); - F_tags.row(fid1) = Vector2i(cur_data, top_data); - } - } - } - } - - // V = V * delta_x; - // igl::writeOFF(output_path, V, F); - // igl::writeOBJ(output_path, V, F); - - logger().info("V = {}, F = {}", V.rows(), F.rows()); -} - -void image_to_tagged_tets(const std::string& filename, MatrixXd& V, MatrixXi& T, VectorXi T_tags) -{ - //// extract surface from image - //{ - // // read raw image data + create triangle soup - // Eigen::MatrixXi F; - // Eigen::MatrixXd V; - // extract_triangle_soup_from_image(filename, F, V); - // igl::writeOFF("triangle_soup_fine.off", V, F); - //} - // - // std::vector verts; - // std::vector> tris; - // std::vector modified_nonmanifold_v; - // std::pair box_minmax; - // - //// read triangle soup - //{ - // const double remove_duplicate_eps = 0.01; - // wmtk::stl_to_manifold_wmtk_input( - // "triangle_soup_fine.off", - // remove_duplicate_eps, - // box_minmax, - // verts, - // tris, - // modified_nonmanifold_v); - //} - // - // std::vector v_simplified; - // std::vector> f_simplified; - // - //// simplification - // app::sec::ShortestEdgeCollapse surf_mesh(verts, 0, false); - //{ - // // must be a small envelope to ensure correct tet tags later on - // surf_mesh.create_mesh(verts.size(), tris, modified_nonmanifold_v, 0.1); - // assert(surf_mesh.check_mesh_connectivity_validity()); - // - // wmtk::logger().info("input {} simplification", filename); - // surf_mesh.collapse_shortest(0); - // surf_mesh.consolidate_mesh(); - // surf_mesh.write_triangle_mesh("triangle_soup_coarse.off"); - // - // //// get the simplified input - // v_simplified.resize(surf_mesh.vert_capacity()); - // f_simplified.resize(surf_mesh.tri_capacity()); - // for (const auto& t : surf_mesh.get_vertices()) { - // const size_t i = t.vid(surf_mesh); - // v_simplified[i] = surf_mesh.vertex_attrs[i].pos; - // } - // - // for (const auto& t : surf_mesh.get_faces()) { - // const auto i = t.fid(surf_mesh); - // const auto vs = surf_mesh.oriented_tri_vids(t); - // for (int j = 0; j < 3; j++) { - // f_simplified[i][j] = vs[j]; - // } - // } - //} - // - // wmtk::remove_duplicates(v_simplified, f_simplified, 0.01); - // - /////////////////////////////////////////////////////// - // - // igl::Timer timer; - // timer.start(); - // - //// triangle insertion with volumeremesher on the simplified mesh - // std::vector v_rational; - // std::vector> facets; - // std::vector is_v_on_input; - // std::vector> tets; - // std::vector tet_face_on_input_surface; - // - // std::cout << "vsimp size: " << v_simplified.size() << std::endl; - // std::cout << "fsimp size: " << f_simplified.size() << std::endl; - // - // igl::Timer insertion_timer; - // insertion_timer.start(); - // - //{ - // ImageSimulationMesh mesh_for_insertion(params, *ptr_env, surf_mesh.m_envelope, - // NUM_THREADS); mesh_for_insertion.insertion_by_volumeremesher( - // v_simplified, - // f_simplified, - // v_rational, - // facets, - // is_v_on_input, - // tets, - // tet_face_on_input_surface); - // } - // - //// generate new mesh - // image_simulation::ImageSimulationMesh mesh(params, *ptr_env, surf_mesh.m_envelope, - // NUM_THREADS); - // - // mesh.init_from_Volumeremesher(v_rational, is_v_on_input, tets, tet_face_on_input_surface); - // - // const double insertion_time = insertion_timer.getElapsedTime(); - // wmtk::logger().info("volume remesher insertion time: {}s", insertion_time); - // - // mesh.consolidate_mesh(); - //{ - // int num_parts = mesh.flood_fill(); - // logger().info("flood fill parts {}", num_parts); - // } - // - //// add tags - // tag_tets_from_image(input_paths[0], mesh); - // - // mesh.write_vtu(params.output_path + "_0.vtu"); -} - -} // namespace wmtk::components::image_simulation \ No newline at end of file diff --git a/components/image_simulation/wmtk/components/image_simulation/extract_soup.hpp b/components/image_simulation/wmtk/components/image_simulation/extract_soup.hpp deleted file mode 100644 index a0b3984468..0000000000 --- a/components/image_simulation/wmtk/components/image_simulation/extract_soup.hpp +++ /dev/null @@ -1,34 +0,0 @@ -#pragma once - -#include -#include - -namespace wmtk::components::image_simulation { - -class ImageSimulationMesh; - -void read_array_data_ascii( - std::vector>>& data, - const std::string& filename); - -void extract_triangle_soup_from_image(std::string filename, Eigen::MatrixXi& F, Eigen::MatrixXd& V); - -/** - * @brief Extracts a triangle soup from a volumetric image. - * - * Whenever two voxels contain different data, two triangles are inserted in between the two voxels. - * - * @param data The voxel image data. - * @param F #Fx3 triangles. - * @param V #Vx3 vertices. - * @param F_tags #Fx2 The voxel data that was separated by this triangle. - */ -void extract_triangle_soup_from_image( - const std::vector>>& data, - MatrixXi& F, - MatrixXd& V, - MatrixXi& F_tags); - -void image_to_tagged_tets(const std::string& filename, MatrixXd& V, MatrixXi& T, VectorXi T_tags); - -} // namespace wmtk::components::image_simulation \ No newline at end of file diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp index 738055091e..59a0afcd90 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp @@ -14,107 +14,32 @@ #include "ImageSimulationMeshTri.hpp" #include "Parameters.h" #include "expression_parser/Parser.hpp" -#include "extract_soup.hpp" #include "read_image_msh.hpp" #include "image_simulation_spec.hpp" namespace wmtk::components::image_simulation { - -void run_3D(const nlohmann::json& json_params, const InputData& input_data) +template +void apply_operation(MeshT& mesh, const nlohmann::json& json_params) { - using wmtk::utils::resolve_path; - using Tuple = TetMesh::Tuple; - - const std::filesystem::path root = - json_params.contains("json_input_file") ? json_params["json_input_file"] : ""; - - Parameters params; - params.output_path = json_params["output"]; - const bool skip_simplify = json_params["skip_simplify"]; - 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"]; - - params.epsr = json_params["eps_rel"]; - params.eps = json_params["eps"]; - params.lr = json_params["length_rel"]; - params.l = json_params["length"]; - params.stop_energy = json_params["stop_energy"]; - params.stop_at_float = json_params["stop_at_float"]; - params.preserve_topology = json_params["preserve_topology"]; - - params.epsr_simplify = json_params["eps_simplify_rel"]; - params.eps_simplify = json_params["eps_simplify"]; - - params.smooth_without_envelope = json_params["smooth_without_envelope"]; - - params.w_amips = json_params["w_amips"]; - - const bool write_vtu = json_params["write_vtu"]; - - params.debug_output = json_params["DEBUG_output"]; - params.perform_sanity_checks = json_params["DEBUG_sanity_checks"]; - - params.operation = json_params["operation"]; - - std::filesystem::path output_filename = params.output_path; - - params.init(input_data.V_input.colwise().minCoeff(), input_data.V_input.colwise().maxCoeff()); - - igl::Timer timer; - timer.start(); - - image_simulation::ImageSimulationMesh mesh(params, params.eps, NUM_THREADS); - // first init envelope - if (input_data.V_envelope.size() != 0) { - mesh.init_envelope(input_data.V_envelope, input_data.F_envelope); + static_assert( + std::is_same_v || + std::is_same_v, + "MeshT must be either ImageSimulationMesh or ImageSimulationMeshTri"); + + if constexpr (std::is_same_v) { + logger().info("Use 3D operation"); + } else if constexpr (std::is_same_v) { + logger().info("Use 2D operation"); } - if (input_data.V_input_r.size() == 0) { - logger().info("Use float input for TetWild"); - mesh.init_from_image( - input_data.V_input, - input_data.T_input, - input_data.T_input_tag, - input_data.tag_names); - } else { - logger().warn("Use RATIONAL input for TetWild"); - mesh.init_from_image( - input_data.V_input_r, - input_data.T_input, - input_data.T_input_tag, - input_data.tag_names); - } - - auto write_unique_vtu = [&write_vtu, &mesh, &output_filename]() { - static size_t vtu_counter = 0; - if (write_vtu) { - mesh.write_vtu(fmt::format("{}_{}", output_filename.string(), vtu_counter++)); - } - }; - - mesh.consolidate_mesh(); - - write_unique_vtu(); - - - if (params.preserve_topology && !skip_simplify) { - // collapse for getting the right edge length - logger().info("Simplify after insertion"); - mesh.simplify(); - } - - const std::string tags_selection_str = json_params["tags_selection"]; - const auto tags_selection_expr = - expression_parser::parse(tags_selection_str, mesh.m_tag_name_to_id); - logger().info("Parsed tags_selection expression: {}", tags_selection_expr->to_string()); - // /////////apply operation - const std::string operation = params.operation; + const std::string operation = json_params["operation"]; if (operation == "remeshing") { - mesh.set_length_regions(json_params["length_region"]); - mesh.mesh_improvement(max_its); // <-- tetwild + if constexpr (std::is_same_v) { + mesh.set_length_regions(json_params["length_region"]); + } + mesh.mesh_improvement(json_params["max_iterations"]); // <-- tetwild } else if (operation == "fill_holes_topo") { const std::vector fill_holes_tags_names = json_params["fill_holes_tags"]; std::vector fill_holes_tags; @@ -218,78 +143,101 @@ void run_3D(const nlohmann::json& json_params, const InputData& input_data) } else { log_and_throw_error("Unknown image simulation operation"); } +} + +void run_3D(const nlohmann::json& json_params, const InputData& input_data) +{ + Parameters params(json_params); + params.init(input_data.V_input.colwise().minCoeff(), input_data.V_input.colwise().maxCoeff()); + + igl::Timer timer; + timer.start(); + + image_simulation::ImageSimulationMesh mesh(params, params.eps, params.NUM_THREADS); + // first init envelope + if (input_data.V_envelope.size() != 0) { + mesh.init_envelope(input_data.V_envelope, input_data.F_envelope); + } + if (input_data.V_input_r.size() == 0) { + logger().info("Use float input for TetWild"); + mesh.init_from_image( + input_data.V_input, + input_data.T_input, + input_data.T_input_tag, + input_data.tag_names); + } else { + logger().warn("Use RATIONAL input for TetWild"); + mesh.init_from_image( + input_data.V_input_r, + input_data.T_input, + input_data.T_input_tag, + input_data.tag_names); + } + + auto write_unique_vtu = [¶ms, &mesh]() { + static size_t vtu_counter = 0; + if (params.write_vtu) { + mesh.write_vtu(fmt::format("{}_{}", params.output_path, vtu_counter++)); + } + }; + + mesh.consolidate_mesh(); + + write_unique_vtu(); + + + if (params.preserve_topology && !params.skip_simplify) { + // collapse for getting the right edge length + logger().info("Simplify after insertion"); + mesh.simplify(); + } + + const std::string tags_selection_str = json_params["tags_selection"]; + const auto tags_selection_expr = + expression_parser::parse(tags_selection_str, mesh.m_tag_name_to_id); + logger().info("Parsed tags_selection expression: {}", tags_selection_expr->to_string()); + + // /////////apply operation + apply_operation(mesh, json_params); mesh.consolidate_mesh(); double time = timer.getElapsedTime(); - wmtk::logger().info("total time {}s", time); + logger().info("total time {}s", time); /////////output auto [max_energy, avg_energy] = mesh.get_max_avg_energy(); - std::string report_filename = json_params["report"]; - if (!report_filename.empty()) { - std::ofstream fout(report_filename); - fout << "#t: " << mesh.tet_size() << std::endl; - fout << "#v: " << mesh.vertex_size() << std::endl; - fout << "max_energy: " << max_energy << std::endl; - fout << "avg_energy: " << avg_energy << std::endl; - fout << "eps: " << params.eps << std::endl; - fout << "threads: " << NUM_THREADS << std::endl; - fout << "time: " << time << std::endl; + const std::string report_file = json_params["report"]; + if (!report_file.empty()) { + std::ofstream fout(report_file); + nlohmann::json report; + report["#t"] = mesh.get_faces().size(); + report["#v"] = mesh.get_vertices().size(); + report["max_energy"] = max_energy; + report["avg_energy"] = avg_energy; + report["eps"] = params.eps; + report["threads"] = params.NUM_THREADS; + report["time"] = time; + fout << std::setw(4) << report; fout.close(); } logger().info("final max energy = {} avg = {}", max_energy, avg_energy); - const bool write_envelope = json_params["write_envelope"]; - // mesh.write_msh(output_filename.string() + ".msh", write_envelope); - mesh.write_msh_groups(output_filename.string() + ".msh", write_envelope); + mesh.write_msh(params.output_path + ".msh", params.write_envelope); write_unique_vtu(); - if (write_vtu) { - mesh.write_surface(output_filename.string() + "_surface.obj"); + if (params.write_vtu) { + mesh.write_surface(params.output_path + "_surface.obj"); } } void run_2D(const nlohmann::json& json_params, const InputData& input_data) { - using wmtk::utils::resolve_path; - using Tuple = TetMesh::Tuple; - - const std::filesystem::path root = - json_params.contains("json_input_file") ? json_params["json_input_file"] : ""; - - Parameters params; - params.output_path = json_params["output"]; - const bool skip_simplify = json_params["skip_simplify"]; - 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"]; - - params.epsr = json_params["eps_rel"]; - params.eps = json_params["eps"]; - params.lr = json_params["length_rel"]; - params.l = json_params["length"]; - params.stop_energy = json_params["stop_energy"]; - params.stop_at_float = json_params["stop_at_float"]; - params.preserve_topology = json_params["preserve_topology"]; - - params.smooth_without_envelope = json_params["smooth_without_envelope"]; - - params.w_amips = json_params["w_amips"]; - - const bool write_vtu = json_params["write_vtu"]; - - params.debug_output = json_params["DEBUG_output"]; - params.perform_sanity_checks = json_params["DEBUG_sanity_checks"]; - - params.operation = json_params["operation"]; - - std::filesystem::path output_filename = params.output_path; - + Parameters params(json_params); params.init(input_data.V_input.colwise().minCoeff(), input_data.V_input.colwise().maxCoeff()); igl::Timer timer; timer.start(); - image_simulation::tri::ImageSimulationMeshTri mesh(params, params.eps, NUM_THREADS); + image_simulation::tri::ImageSimulationMeshTri mesh(params, params.eps, params.NUM_THREADS); // first init envelope if (input_data.V_envelope.size() != 0) { mesh.init_envelope(input_data.V_envelope, input_data.F_envelope); @@ -303,155 +251,48 @@ void run_2D(const nlohmann::json& json_params, const InputData& input_data) input_data.T_input_tag, input_data.tag_names); - auto write_unique_vtu = [&write_vtu, &mesh, &output_filename]() { + auto write_unique_vtu = [¶ms, &mesh]() { static size_t vtu_counter = 0; - if (write_vtu) { - mesh.write_vtu(fmt::format("{}_{}", output_filename.string(), vtu_counter++)); + if (params.write_vtu) { + mesh.write_vtu(fmt::format("{}_{}", params.output_path, vtu_counter++)); } }; mesh.consolidate_mesh(); - write_unique_vtu(); // /////////apply operation - const std::string operation = params.operation; - if (operation == "remeshing") { - mesh.mesh_improvement(max_its); // <-- tetwild - } else if (operation == "fill_holes_topo") { - const std::vector fill_holes_tags_names = json_params["fill_holes_tags"]; - std::vector fill_holes_tags; - for (const auto& tag_set_names : fill_holes_tags_names) { - // fill_holes_tags.push_back(mesh.string_set_to_cell_tag(tag_set_names)); - const auto expr = expression_parser::parse(tag_set_names, mesh.m_tag_name_to_id); - logger().info("Parsed fill_holes_tags expression: {}", expr->to_string()); - if (!expr->contains_only_and()) { - log_and_throw_error("Only AND operation is allowed in fill_holes_tags expression."); - } - fill_holes_tags.push_back(expr->tags_involved()); - } - const double raw_threshold = json_params["fill_holes_threshold"]; - const double threshold = - raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; - mesh.fill_holes_topo(fill_holes_tags, threshold); - } else if (operation == "tight_seal_topo") { - // tight_seal_tag_sets is a list of lists: [[t1,t2],[t3,...]] - const std::vector> tag_sets_names = - json_params["tight_seal_tag_sets"]; - std::vector> tag_sets; - for (const auto& tag_set_names_vec : tag_sets_names) { - std::vector tag_set_vec; - for (const auto& tag_set_names : tag_set_names_vec) { - auto expr = expression_parser::parse(tag_set_names, mesh.m_tag_name_to_id); - logger().info("Parsed tight_seal_tag_sets expression: {}", expr->to_string()); - if (!expr->contains_only_and()) { - log_and_throw_error( - "Only AND operation is allowed in tight_seal_tag_sets expression."); - } - auto tag_set = expr->tags_involved(); - tag_set_vec.push_back(tag_set); - } - tag_sets.push_back(tag_set_vec); - } - const double raw_threshold = json_params["tight_seal_threshold"]; - const double threshold = - raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; - mesh.tight_seal_topo(tag_sets, threshold); - } else if (operation == "keep_lcc") { - const std::vector lcc_tags_names = json_params["keep_lcc_tags"]; - std::vector lcc_tags; - for (const auto& tag_set_names : lcc_tags_names) { - auto expr = expression_parser::parse(tag_set_names, mesh.m_tag_name_to_id); - logger().info("Parsed keep_lcc_tags expression: {}", expr->to_string()); - if (!expr->contains_only_and()) { - log_and_throw_error("Only AND operation is allowed in keep_lcc_tags expression."); - } - lcc_tags.push_back(expr->tags_involved()); - } - const size_t n_lcc = json_params["keep_lcc_num"]; - mesh.keep_largest_connected_component(lcc_tags, n_lcc); - } else if (operation == "resolve_intersections") { - const std::vector tags_names = json_params["resolve_intersections_tags"]; - std::vector tags; - for (const auto& tag_set_names : tags_names) { - auto expr = expression_parser::parse(tag_set_names, mesh.m_tag_name_to_id); - logger().info("Parsed resolve_intersections_tags expression: {}", expr->to_string()); - if (!expr->contains_only_and()) { - log_and_throw_error( - "Only AND operation is allowed in resolve_intersections_tags expression."); - } - auto tag_set = expr->tags_involved(); - if (tag_set.size() != 2) { - log_and_throw_error( - "Exactly two tags must be provided in each resolve_intersections_tags " - "expression. Expression '{}' involves {} tags.", - expr->to_string(), - tag_set.size()); - } - tags.push_back(tag_set); - } - mesh.resolve_intersections(tags); - } else if (operation == "replace_tags") { - const std::vector tags_in_names = json_params["replace_tags_in"]; - std::vector tags_in; - for (const auto& tag_set_names : tags_in_names) { - const auto expr_in = expression_parser::parse(tag_set_names, mesh.m_tag_name_to_id); - logger().info("Parsed tags_in expression: {}", expr_in->to_string()); - if (!expr_in->contains_only_and()) { - log_and_throw_error("Only AND operation is allowed in replace_tags_in expression."); - } - tags_in.push_back(expr_in->tags_involved()); - } - const std::string tag_out_name = json_params["replace_tags_out"]; - const auto expr_out = expression_parser::parse(tag_out_name, mesh.m_tag_name_to_id); - logger().info("Parsed replace_tags_out expression: {}", expr_out->to_string()); - if (!expr_out->contains_only_and()) { - log_and_throw_error("Only AND operation is allowed in replace_tags_out expression."); - } - CellTag tag_out = mesh.string_set_to_cell_tag(expr_out->tag_names_involved()); - mesh.replace_tags(tags_in, tag_out); - } else if (operation == "tag_priority") { - const std::vector tag_priority_names = json_params["tag_priority"]; - std::vector tag_priority; - for (const auto& tag_name : tag_priority_names) { - CellTag tag = mesh.string_set_to_cell_tag({tag_name}); - tag_priority.push_back(*tag.begin()); - } - mesh.tag_priority(tag_priority); - } else { - log_and_throw_error("Unknown image simulation operation"); - } + apply_operation(mesh, json_params); mesh.consolidate_mesh(); double time = timer.getElapsedTime(); - wmtk::logger().info("total time {}s", time); + logger().info("total time {}s", time); /////////output auto [max_energy, avg_energy] = mesh.get_max_avg_energy(); - std::string report_filename = json_params["report"]; - if (!report_filename.empty()) { - std::ofstream fout(report_filename); - fout << "#t: " << mesh.get_faces().size() << std::endl; - fout << "#v: " << mesh.get_vertices().size() << std::endl; - fout << "max_energy: " << max_energy << std::endl; - fout << "avg_energy: " << avg_energy << std::endl; - fout << "eps: " << params.eps << std::endl; - fout << "threads: " << NUM_THREADS << std::endl; - fout << "time: " << time << std::endl; + const std::string report_file = json_params["report"]; + if (!report_file.empty()) { + std::ofstream fout(report_file); + nlohmann::json report; + report["#t"] = mesh.get_faces().size(); + report["#v"] = mesh.get_vertices().size(); + report["max_energy"] = max_energy; + report["avg_energy"] = avg_energy; + report["eps"] = params.eps; + report["threads"] = params.NUM_THREADS; + report["time"] = time; + fout << std::setw(4) << report; fout.close(); } - wmtk::logger().info("final max energy = {} avg = {}", max_energy, avg_energy); - const bool write_envelope = json_params["write_envelope"]; - // mesh.write_msh(output_filename.string() + ".msh", write_envelope); - mesh.write_msh_groups(output_filename.string() + ".msh", write_envelope); + logger().info("final max energy = {} avg = {}", max_energy, avg_energy); + mesh.write_msh(params.output_path + ".msh", params.write_envelope); write_unique_vtu(); } void image_simulation(nlohmann::json json_params) { using wmtk::utils::resolve_path; - using Tuple = TetMesh::Tuple; // verify input and inject defaults { @@ -502,8 +343,6 @@ void image_simulation(nlohmann::json json_params) std::string extension = std::filesystem::path(input_paths[0]).extension().string(); if (extension == ".msh") { input_data = read_image_msh(input_paths[0]); - } else if (extension == ".raw") { - input_data = read_image(input_paths, output_filename.string(), json_params); } else { input_data = read_mesh(input_paths, output_filename.string(), json_params); } @@ -513,7 +352,7 @@ void image_simulation(nlohmann::json json_params) } else { run_2D(json_params, input_data); } - wmtk::logger().info("======= finish ========="); + logger().info("======= finish ========="); } } // namespace wmtk::components::image_simulation \ No newline at end of file diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp index d90f8a3307..763157297d 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp @@ -12,11 +12,9 @@ nlohmann::json image_simulation_spec = R"( "operation", "output", "input_names", - "ijk_to_ras", "input_transform", "skip_simplify", "use_sample_envelope", - "use_tetgen", "num_threads", "max_iterations", "eps_rel", @@ -30,7 +28,6 @@ nlohmann::json image_simulation_spec = R"( "stop_at_float", "preserve_topology", "w_amips", - "smooth_without_envelope", "write_vtu", "write_envelope", "log_file", @@ -99,29 +96,6 @@ nlohmann::json image_simulation_spec = R"( "default": "out", "doc": "Output file name (without extension)." }, - { - "pointer": "/ijk_to_ras", - "type": "list", - "doc": "Transformation matrix (4x4 homogeneous coordinates) from image coordinates to the RAS coordinate system.", - "min": 4, - "max": 4, - "default": [ - [1, 0, 0, 0], - [0, 1, 0, 0], - [0, 0, 1, 0], - [0, 0, 0, 1] - ] - }, - { - "pointer": "/ijk_to_ras/*", - "type": "list", - "min": 4, - "max": 4 - }, - { - "pointer": "/ijk_to_ras/*/*", - "type": "float" - }, { "pointer": "/input_transform", "type": "list", @@ -156,12 +130,6 @@ nlohmann::json image_simulation_spec = R"( "default": false, "doc": "Use sample envelope instead of exact one." }, - { - "pointer": "/use_tetgen", - "type": "bool", - "default": false, - "doc": "Use tetgen for insertion instead of binary space partitioning." - }, { "pointer": "/num_threads", "type": "int", @@ -266,12 +234,6 @@ nlohmann::json image_simulation_spec = R"( "max": 1, "doc": "AMIPS energy" }, - { - "pointer": "/smooth_without_envelope", - "type": "bool", - "default": false, - "doc": "Ignore the envelope during smoothing and re-compute it afterwards. This may cause the geometry to drift!" - }, { "pointer": "/write_vtu", "type": "bool", diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json index 523d7861c6..e5162d605b 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json @@ -7,11 +7,9 @@ "operation", "output", "input_names", - "ijk_to_ras", "input_transform", "skip_simplify", "use_sample_envelope", - "use_tetgen", "num_threads", "max_iterations", "eps_rel", @@ -25,7 +23,6 @@ "stop_at_float", "preserve_topology", "w_amips", - "smooth_without_envelope", "write_vtu", "write_envelope", "log_file", @@ -94,29 +91,6 @@ "default": "out", "doc": "Output file name (without extension)." }, - { - "pointer": "/ijk_to_ras", - "type": "list", - "doc": "Transformation matrix (4x4 homogeneous coordinates) from image coordinates to the RAS coordinate system.", - "min": 4, - "max": 4, - "default": [ - [1, 0, 0, 0], - [0, 1, 0, 0], - [0, 0, 1, 0], - [0, 0, 0, 1] - ] - }, - { - "pointer": "/ijk_to_ras/*", - "type": "list", - "min": 4, - "max": 4 - }, - { - "pointer": "/ijk_to_ras/*/*", - "type": "float" - }, { "pointer": "/input_transform", "type": "list", @@ -151,12 +125,6 @@ "default": false, "doc": "Use sample envelope instead of exact one." }, - { - "pointer": "/use_tetgen", - "type": "bool", - "default": false, - "doc": "Use tetgen for insertion instead of binary space partitioning." - }, { "pointer": "/num_threads", "type": "int", @@ -261,12 +229,6 @@ "max": 1, "doc": "AMIPS energy" }, - { - "pointer": "/smooth_without_envelope", - "type": "bool", - "default": false, - "doc": "Ignore the envelope during smoothing and re-compute it afterwards. This may cause the geometry to drift!" - }, { "pointer": "/write_vtu", "type": "bool", diff --git a/components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp b/components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp index c0cb29058a..ff5d50d514 100644 --- a/components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp @@ -10,32 +10,6 @@ namespace wmtk::components::image_simulation { -Matrix4d get_ijk2xyz(const nlohmann::json& json_params) -{ - const std::array, 4> ijk_to_ras = json_params["ijk_to_ras"]; - Matrix4d ijk2ras; - for (size_t i = 0; i < 4; ++i) { - for (size_t j = 0; j < 4; ++j) { - ijk2ras(i, j) = ijk_to_ras[i][j]; - } - } - logger().info("IJK to RAS:\n{}", ijk2ras); - - /** - * R = -x - * A = z - * S = y - */ - Matrix4d ras2xyz; - ras2xyz << -1, 0, 0, 0, // - 0, 0, 1, 0, // - 0, 1, 0, 0, // - 0, 0, 0, 1; - Matrix4d ijk2xyz = ras2xyz * ijk2ras; - - return ijk2xyz; -} - Matrix4d get_input_transform(const nlohmann::json& json_params, const size_t input_index) { const nlohmann::json& its_j = json_params["input_transform"]; @@ -403,63 +377,6 @@ InputData read_image_msh(const std::string& path) return input_data; } -InputData read_image( - const std::vector& input_paths, - const std::string& output_filename, - const nlohmann::json& json_params) -{ - InputData input_data; - - Matrix4d ijk2xyz = get_ijk2xyz(json_params); - - 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 Vector4d single_voxel_max = ijk2xyz * Vector4d::Ones(); - const Vector4d single_voxel_min = ijk2xyz * Vector4d(0, 0, 0, 1); - double eps = (from_homogenuous(single_voxel_max) - from_homogenuous(single_voxel_min)) - .cwiseAbs() - .minCoeff() * - 0.1; - if (eps <= 0) { - logger().warn("EPS = {}, ijk_to_ras matix might be broken! Changing eps to 1e-4", eps); - eps = 1e-4; - } - - // convert image into tet mesh - EmbedSurface image_mesh(input_paths, ijk2xyz); - - if (write_vtu) { - image_mesh.write_surf_off(output_filename + "_input.off"); - } - - // input_data.V_envelope = image_mesh.V_surface(); - // input_data.F_envelope = image_mesh.F_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) { - // image_mesh.write_emb_vtu(get_unique_vtu_name()); - image_mesh.write_emb_surf_off(output_filename + "_input_emb.off"); - } - - input_data.V_input = image_mesh.V_emb(); - if (!all_rounded) { - input_data.V_input_r = image_mesh.V_emb_r(); - } - input_data.T_input = image_mesh.T_emb(); - input_data.T_input_tag = image_mesh.T_tags(); - - wmtk::logger().info("======= finish image-tet conversion ========="); - - return input_data; -} - InputData read_mesh( const std::vector& input_paths, const std::string& output_filename, @@ -519,9 +436,7 @@ InputData read_mesh( input_data.V_envelope = image_mesh.V_surface(); input_data.F_envelope = image_mesh.F_surface(); - const bool all_rounded = json_params["use_tetgen"] - ? image_mesh.embed_surface_tetgen() - : image_mesh.embed_surface(preserve_topology); + const bool all_rounded = image_mesh.embed_surface(preserve_topology); image_mesh.consolidate(); if (debug_output) { diff --git a/components/image_simulation/wmtk/components/image_simulation/read_image_msh.hpp b/components/image_simulation/wmtk/components/image_simulation/read_image_msh.hpp index 27981bd8f7..c54cf0bd0f 100644 --- a/components/image_simulation/wmtk/components/image_simulation/read_image_msh.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/read_image_msh.hpp @@ -23,16 +23,6 @@ struct InputData */ InputData read_image_msh(const std::string& path); -/** - * @brief Read one or multiple images and convert them into a tet mesh. - * - * The color in the image become tet tags. - */ -InputData read_image( - const std::vector& input_paths, - const std::string& output_filename, - const nlohmann::json& json_params); - /** * @brief Read one or multiple meshes and convert them into a tet mesh. * From 4fd30b320fcf6a13c7d4b18a9e60ede4b4f2ed8e Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Thu, 11 Jun 2026 15:12:10 +0200 Subject: [PATCH 2/3] replace_tags has a list of output tags. This enables the replacement of multiple tags in one run. --- cmake/wmtk_data.cmake | 2 +- .../image_simulation/AnnotationsTet.cpp | 29 ++++++++++++------- .../image_simulation/AnnotationsTri.cpp | 27 ++++++++++------- .../image_simulation/EmbedSurface.cpp | 8 ++--- .../image_simulation/EmbedSurface.hpp | 1 + .../image_simulation/ImageSimulationMesh.h | 2 +- .../ImageSimulationMeshTri.hpp | 2 +- .../image_simulation/image_simulation.cpp | 25 +++++++++++----- .../image_simulation_spec.hpp | 9 ++++-- .../image_simulation_spec.json | 9 ++++-- .../image_simulation/read_image_msh.cpp | 9 +++--- 11 files changed, 78 insertions(+), 45 deletions(-) diff --git a/cmake/wmtk_data.cmake b/cmake/wmtk_data.cmake index fbfa14508e..33c2b4e2f2 100644 --- a/cmake/wmtk_data.cmake +++ b/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 a97b2a974e80cc288624f86c9cc03eae7f050c4f + GIT_TAG f2bd8ed5e8af5a125ee4d92cd0af8f32cec60e02 CONFIGURE_COMMAND "" BUILD_COMMAND "" diff --git a/components/image_simulation/wmtk/components/image_simulation/AnnotationsTet.cpp b/components/image_simulation/wmtk/components/image_simulation/AnnotationsTet.cpp index 4d4166a281..631164790c 100644 --- a/components/image_simulation/wmtk/components/image_simulation/AnnotationsTet.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/AnnotationsTet.cpp @@ -579,22 +579,31 @@ void ImageSimulationMesh::resolve_intersections(const std::vector& inte m_envelope.reset(); } -void ImageSimulationMesh::replace_tags(const std::vector& tags_in, const CellTag& tag_out) +void ImageSimulationMesh::replace_tags( + const std::vector& tags_in, + const std::vector& tags_out) { for (const Tuple& t : get_tets()) { CellTag& tags = m_tet_attribute[t.tid(*this)].tags; - bool found_some = false; - for (const CellTag& tag : tags_in) { - if (set_includes(tags, tag)) { - found_some = true; - // remove tag from tags - for (const size_t t : tag) { - tags.erase(t); + std::vector found_tags(tags_in.size(), false); + for (size_t i = 0; i < tags_in.size(); ++i) { + if (set_includes(tags, tags_in[i])) { + found_tags[i] = true; + } + } + // remove "tags_in" from tags + for (size_t i = 0; i < tags_in.size(); ++i) { + if (found_tags[i]) { + for (const int64_t tt : tags_in[i]) { + tags.erase(tt); } } } - if (found_some) { - tags.insert(tag_out.begin(), tag_out.end()); + // add "tags_out" to tags + for (size_t i = 0; i < tags_out.size(); ++i) { + if (found_tags[i]) { + tags.insert(tags_out[i].begin(), tags_out[i].end()); + } } } diff --git a/components/image_simulation/wmtk/components/image_simulation/AnnotationsTri.cpp b/components/image_simulation/wmtk/components/image_simulation/AnnotationsTri.cpp index b01bc3adc1..6dc268482b 100644 --- a/components/image_simulation/wmtk/components/image_simulation/AnnotationsTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/AnnotationsTri.cpp @@ -535,22 +535,29 @@ void ImageSimulationMeshTri::resolve_intersections(const std::vector& i void ImageSimulationMeshTri::replace_tags( const std::vector& tags_in, - const CellTag& tag_out) + const std::vector& tags_out) { for (const Tuple& t : get_faces()) { CellTag& tags = m_face_attribute[t.fid(*this)].tags; - bool found_some = false; - for (const CellTag& tag : tags_in) { - if (set_includes(tags, tag)) { - found_some = true; - // remove tag from tags - for (const size_t t : tag) { - tags.erase(t); + std::vector found_tags(tags_in.size(), false); + for (size_t i = 0; i < tags_in.size(); ++i) { + if (set_includes(tags, tags_in[i])) { + found_tags[i] = true; + } + } + // remove "tags_in" from tags + for (size_t i = 0; i < tags_in.size(); ++i) { + if (found_tags[i]) { + for (const int64_t tt : tags_in[i]) { + tags.erase(tt); } } } - if (found_some) { - tags.insert(tag_out.begin(), tag_out.end()); + // add "tags_out" to tags + for (size_t i = 0; i < tags_out.size(); ++i) { + if (found_tags[i]) { + tags.insert(tags_out[i].begin(), tags_out[i].end()); + } } } diff --git a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp index 2ddce88bc6..9acf86f555 100644 --- a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp @@ -1,11 +1,7 @@ #include "EmbedSurface.hpp" // clang-format off -#include -#include -#include #include -#include #include #include #include @@ -792,7 +788,7 @@ bool EmbedSurface::embed_surface(const bool flood_fill) } // get all tags for one image - for (size_t img_id = 0; img_id < m_T_tags.cols(); ++img_id) { + for (size_t img_id = 0; img_id < Fs.size(); ++img_id) { std::map m; for (size_t j = 0; j < region.size(); ++j) { int tag = m_T_tags.coeff(region[j], img_id); @@ -984,6 +980,7 @@ void EmbedSurface::tag_from_winding_number() { m_T_tags.resize(m_T_emb.rows(), Fs.size()); m_T_tags.setZero(); + m_tags.resize(m_T_emb.rows()); MatrixXd P; P.resize(m_T_emb.rows(), 3); @@ -1005,6 +1002,7 @@ void EmbedSurface::tag_from_winding_number() continue; } m_T_tags.coeffRef(j, i) = 1; + m_tags[j].insert(i); } } m_T_tags.makeCompressed(); diff --git a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp index c437bc3d4f..61cf51ed6d 100644 --- a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp @@ -141,6 +141,7 @@ class EmbedSurface MatrixXi m_F_on_surface; // tags on the tets MatrixSi m_T_tags; + std::vector> m_tags; // input from triangle meshes std::vector Vs; diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h index 979dae0bae..bd61af28f4 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h @@ -608,7 +608,7 @@ class ImageSimulationMesh : public wmtk::TetMesh void resolve_intersections(const std::vector& intersecting_tags); - void replace_tags(const std::vector& tags_in, const CellTag& tag_out); + void replace_tags(const std::vector& tags_in, const std::vector& tags_out); /** * @brief Gives tags priority over others. diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index 39f4e128c2..99c0751408 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp @@ -311,7 +311,7 @@ class ImageSimulationMeshTri : public wmtk::TriMesh void resolve_intersections(const std::vector& intersecting_tags); - void replace_tags(const std::vector& tags_in, const CellTag& tag_out); + void replace_tags(const std::vector& tags_in, const std::vector& tags_out); void tag_priority(const std::vector& tags_order); diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp index 59a0afcd90..cbe3ebfe48 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp @@ -115,6 +115,11 @@ void apply_operation(MeshT& mesh, const nlohmann::json& json_params) mesh.resolve_intersections(tags); } else if (operation == "replace_tags") { const std::vector tags_in_names = json_params["replace_tags_in"]; + const std::vector tag_out_names = json_params["replace_tags_out"]; + if (tags_in_names.size() != tag_out_names.size()) { + log_and_throw_error( + "replace_tags_in and replace_tags_out must have the same number of expressions."); + } std::vector tags_in; for (const auto& tag_set_names : tags_in_names) { const auto expr_in = expression_parser::parse(tag_set_names, mesh.m_tag_name_to_id); @@ -124,14 +129,18 @@ void apply_operation(MeshT& mesh, const nlohmann::json& json_params) } tags_in.push_back(expr_in->tags_involved()); } - const std::string tag_out_name = json_params["replace_tags_out"]; - const auto expr_out = expression_parser::parse(tag_out_name, mesh.m_tag_name_to_id); - logger().info("Parsed replace_tags_out expression: {}", expr_out->to_string()); - if (!expr_out->contains_only_and()) { - log_and_throw_error("Only AND operation is allowed in replace_tags_out expression."); + std::vector tags_out; + for (const auto& tag_out_name : tag_out_names) { + const auto expr_out = expression_parser::parse(tag_out_name, mesh.m_tag_name_to_id); + logger().info("Parsed replace_tags_out expression: {}", expr_out->to_string()); + if (!expr_out->contains_only_and()) { + log_and_throw_error( + "Only AND operation is allowed in replace_tags_out expression."); + } + CellTag tag_out = mesh.string_set_to_cell_tag(expr_out->tag_names_involved()); + tags_out.push_back(tag_out); } - CellTag tag_out = mesh.string_set_to_cell_tag(expr_out->tag_names_involved()); - mesh.replace_tags(tags_in, tag_out); + mesh.replace_tags(tags_in, tags_out); } else if (operation == "tag_priority") { const std::vector tag_priority_names = json_params["tag_priority"]; std::vector tag_priority; @@ -186,7 +195,7 @@ void run_3D(const nlohmann::json& json_params, const InputData& input_data) write_unique_vtu(); - if (params.preserve_topology && !params.skip_simplify) { + if (params.operation == "remeshing" && params.preserve_topology && !params.skip_simplify) { // collapse for getting the right edge length logger().info("Simplify after insertion"); mesh.simplify(); diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp index 763157297d..24445da023 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp @@ -352,9 +352,14 @@ nlohmann::json image_simulation_spec = R"( }, { "pointer": "/replace_tags_out", + "type": "list", + "default": [], + "doc": "The output tags (written as intersection) that should replace the input tags. For each input tag, an output tag must be provided. The order of the output tags should correspond to the order of the input tags." + }, + { + "pointer": "/replace_tags_out/*", "type": "string", - "default": "_", - "doc": "The output tags (written as intersection) that should replace the input tags." + "doc": "An intersection of tags that should replace the input tags." }, { "pointer": "/tag_priority", diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json index e5162d605b..5fe8d9f658 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json @@ -347,9 +347,14 @@ }, { "pointer": "/replace_tags_out", + "type": "list", + "default": [], + "doc": "The output tags (written as intersection) that should replace the input tags. For each input tag, an output tag must be provided. The order of the output tags should correspond to the order of the input tags." + }, + { + "pointer": "/replace_tags_out/*", "type": "string", - "default": "_", - "doc": "The output tags (written as intersection) that should replace the input tags." + "doc": "An intersection of tags that should replace the input tags." }, { "pointer": "/tag_priority", diff --git a/components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp b/components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp index ff5d50d514..c794e07c6c 100644 --- a/components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp @@ -389,9 +389,7 @@ InputData read_mesh( input_transforms[i] = get_input_transform(json_params, i); } - 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 debug_output = json_params["DEBUG_output"]; const bool preserve_topology = json_params["preserve_topology"]; const bool skip_simplify = json_params["skip_simplify"]; @@ -416,8 +414,8 @@ InputData read_mesh( 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. + * 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(); @@ -449,10 +447,11 @@ InputData read_mesh( } input_data.T_input = image_mesh.T_emb(); input_data.T_input_tag = image_mesh.T_tags(); + // create generic tag names first for (int i = 0; i < input_data.T_input_tag.cols(); ++i) { input_data.tag_names.push_back("tag_" + std::to_string(i)); } - + // use provided tag names if possible for (int i = 0; i < input_data.T_input_tag.cols(); ++i) { if (i < input_names.size()) { input_data.tag_names[i] = input_names[i]; From d4954b671cd7034161147b409ce2d05b8a75ccb8 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Fri, 12 Jun 2026 09:27:15 +0200 Subject: [PATCH 3/3] Add missing include. --- .../wmtk/components/image_simulation/EmbedSurface.cpp | 2 +- .../wmtk/components/image_simulation/EmbedSurface.hpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp index 9acf86f555..d984c0aa8b 100644 --- a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp @@ -1002,7 +1002,7 @@ void EmbedSurface::tag_from_winding_number() continue; } m_T_tags.coeffRef(j, i) = 1; - m_tags[j].insert(i); + m_tags[j].insert((int64_t)i); } } m_T_tags.makeCompressed(); diff --git a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp index 61cf51ed6d..112cedc4a1 100644 --- a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include