Skip to content

Commit e7b8bf1

Browse files
daniel-zintdaniel-zint
committed
Add tetgen call to image simulation.
1 parent 9cf65c1 commit e7b8bf1

7 files changed

Lines changed: 141 additions & 11 deletions

File tree

‎cmake/recipes/libigl.cmake‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ CPMAddPackage(
1212
GIT_TAG 3ea7f9480967fcf6bf02ce9b993c0ea6d2fc45f6
1313
OPTIONS
1414
LIBIGL_PREDICATES ON
15+
LIBIGL_COPYLEFT_TETGEN ON
1516
)
1617
# include(eigen)
1718
FetchContent_MakeAvailable(libigl)

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,7 @@ target_link_libraries(wmtk_${COMPONENT_NAME} PUBLIC
4545
wmtk::sec_lib
4646
VolumeRemesher::VolumeRemesher
4747
jse::jse
48+
igl_copyleft::tetgen
4849
)
4950

5051
######################

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

Lines changed: 120 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include <igl/is_vertex_manifold.h>
66
#include <igl/remove_duplicate_vertices.h>
77
#include <igl/remove_unreferenced.h>
8+
#include <igl/copyleft/tetgen/tetrahedralize.h>
89
// igl must be included BEFORE VolumeRemesher
910
#include <VolumeRemesher/embed.h>
1011
#include <VolumeRemesher/numerics.h>
@@ -784,6 +785,8 @@ void EmbedSurface::remove_duplicates(const double eps)
784785

785786
bool EmbedSurface::embed_surface()
786787
{
788+
logger().info("Embed with VolumeInsertion");
789+
787790
std::shared_ptr<SampleEnvelope> ptr_env;
788791
{
789792
const auto v_simplified = V_surf_to_vector();
@@ -838,6 +841,123 @@ bool EmbedSurface::embed_surface()
838841
return all_rounded;
839842
}
840843

844+
bool EmbedSurface::embed_surface_tetgen()
845+
{
846+
logger().info("Embed with TetGen");
847+
848+
std::shared_ptr<SampleEnvelope> ptr_env;
849+
{
850+
const auto v_simplified = V_surf_to_vector();
851+
852+
std::vector<Eigen::Vector3i> tempF(m_F_surface.rows());
853+
for (size_t i = 0; i < tempF.size(); ++i) {
854+
tempF[i] = m_F_surface.row(i);
855+
}
856+
ptr_env = std::make_shared<SampleEnvelope>();
857+
ptr_env->init(v_simplified, tempF, 0.5);
858+
}
859+
860+
Eigen::MatrixXd V_in;
861+
Eigen::MatrixXi F_in;
862+
{
863+
const Vector3d vertices_max = m_V_surface.colwise().maxCoeff();
864+
const Vector3d vertices_min = m_V_surface.colwise().minCoeff();
865+
const double diag = (vertices_max - vertices_min).norm();
866+
867+
// bbox
868+
double delta = diag / 15.0;
869+
const Vector3d box_min(
870+
vertices_min[0] - delta,
871+
vertices_min[1] - delta,
872+
vertices_min[2] - delta);
873+
const Vector3d box_max(
874+
vertices_max[0] + delta,
875+
vertices_max[1] + delta,
876+
vertices_max[2] + delta);
877+
878+
879+
// add corners of domain
880+
std::vector<Vector3d> points(8);
881+
for (int i = 0; i < 8; i++) {
882+
Vector3d& p = points[i];
883+
std::bitset<sizeof(int) * 8> a(i);
884+
for (int j = 0; j < 3; j++) {
885+
if (a.test(j)) {
886+
p[j] = box_max[j];
887+
} else {
888+
p[j] = box_min[j];
889+
}
890+
}
891+
}
892+
893+
//const double voxel_resolution = diag / 20.0;
894+
//std::array<int, 3> N; // number of grid points per dimension
895+
//std::array<double, 3> h; // distance between grid points per dimension
896+
// for (int i = 0; i < 3; i++) {
897+
// const double D = box_max[i] - box_min[i];
898+
// N[i] = (D / voxel_resolution) + 1;
899+
// h[i] = D / N[i];
900+
//}
901+
//
902+
// std::array<std::vector<double>, 3> ds;
903+
// for (int i = 0; i < 3; i++) {
904+
// ds[i].push_back(box_min[i]);
905+
// for (int j = 0; j < N[i] - 1; j++) {
906+
// ds[i].push_back(box_min[i] + h[i] * (j + 1));
907+
// }
908+
// ds[i].push_back(box_max[i]);
909+
//}
910+
//
911+
//const double min_dis = voxel_resolution * voxel_resolution / 4;
912+
//// double min_dis = state.target_edge_len * state.target_edge_len;//epsilon*2
913+
// for (int i = 0; i < ds[0].size(); i++) {
914+
// for (int j = 0; j < ds[1].size(); j++) {
915+
// for (int k = 0; k < ds[2].size(); k++) {
916+
// if ((i == 0 || i == ds[0].size() - 1) && (j == 0 || j == ds[1].size() - 1) &&
917+
// (k == 0 || k == ds[2].size() - 1)) {
918+
// continue;
919+
// }
920+
// const Vector3d p(ds[0][i], ds[1][j], ds[2][k]);
921+
//
922+
// Eigen::Vector3d n;
923+
// const double sqd = ptr_env->nearest_point(p, n);
924+
//
925+
// if (sqd < min_dis) {
926+
// continue;
927+
// }
928+
// points.emplace_back(ds[0][i], ds[1][j], ds[2][k]);
929+
// }
930+
// }
931+
// }
932+
933+
934+
V_in.resize(m_V_surface.rows() + points.size(), 3);
935+
V_in.block(0, 0, m_V_surface.rows(), 3) = m_V_surface;
936+
for (size_t i = 0; i < points.size(); ++i) {
937+
V_in.row(m_V_surface.rows() + i) = points[i];
938+
}
939+
940+
F_in = m_F_surface;
941+
}
942+
943+
MatrixXi TF;
944+
945+
int ret = igl::copyleft::tetgen::tetrahedralize(V_in, F_in, "pq1.414Yc", m_V_emb, m_T_emb, TF);
946+
if (ret != 0) {
947+
log_and_throw_error("Tetwild returned with {}", ret);
948+
}
949+
950+
m_V_emb_r.resizeLike(m_V_emb);
951+
for (int i = 0; i < m_V_emb.rows(); ++i) {
952+
m_V_emb_r.row(i) = to_rational(m_V_emb.row(i));
953+
}
954+
955+
// add tags
956+
tag_tets_from_images(m_img_datas, m_xyz2ijk, m_V_emb, m_T_emb, m_T_tags);
957+
958+
return true;
959+
}
960+
841961
void EmbedSurface::consolidate()
842962
{
843963
std::map<size_t, size_t> old2new;

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

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -110,6 +110,8 @@ class EmbedSurface
110110

111111
bool embed_surface();
112112

113+
bool embed_surface_tetgen();
114+
113115
/**
114116
* @brief Remove unreferenced vertices.
115117
*/

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

Lines changed: 3 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -6,20 +6,10 @@
66

77
#include <jse/jse.h>
88
#include <wmtk/TetMesh.h>
9-
#include <wmtk/utils/Partitioning.h>
10-
#include <paraviewo/VTUWriter.hpp>
119
#include <wmtk/Types.hpp>
12-
#include <wmtk/envelope/Envelope.hpp>
13-
#include <wmtk/utils/InsertTriangleUtils.hpp>
1410
#include <wmtk/utils/Logger.hpp>
15-
#include <wmtk/utils/ManifoldUtils.hpp>
16-
#include <wmtk/utils/Reader.hpp>
17-
#include <wmtk/utils/io.hpp>
18-
#include <wmtk/utils/partition_utils.hpp>
1911
#include <wmtk/utils/resolve_path.hpp>
2012

21-
#include <sec/ShortestEdgeCollapse.h>
22-
2313
#include "EmbedSurface.hpp"
2414
#include "ImageSimulationMesh.h"
2515
#include "Parameters.h"
@@ -164,7 +154,9 @@ void image_simulation(nlohmann::json json_params)
164154
V_envelope = image_mesh.V_surface();
165155
F_envelope = image_mesh.F_surface();
166156

167-
const bool all_rounded = image_mesh.embed_surface();
157+
// const bool all_rounded = image_mesh.embed_surface();
158+
const bool all_rounded = json_params["use_tetgen"] ? image_mesh.embed_surface_tetgen()
159+
: image_mesh.embed_surface();
168160
image_mesh.consolidate();
169161

170162
if (write_vtu) {

‎components/tet_remeshing/wmtk/components/tet_remeshing/tet_remeshing_spec.hpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@ nlohmann::json tet_remeshing_spec = R"(
1313
"ijk_to_ras",
1414
"skip_simplify",
1515
"use_sample_envelope",
16+
"use_tetgen",
1617
"num_threads",
1718
"max_iterations",
1819
"eps_rel",
@@ -82,6 +83,12 @@ nlohmann::json tet_remeshing_spec = R"(
8283
"default": false,
8384
"doc": "Use sample envelope instead of exact one."
8485
},
86+
{
87+
"pointer": "/use_tetgen",
88+
"type": "bool",
89+
"default": false,
90+
"doc": "Use tetgen for embedding an image. Potentially faster but could fail."
91+
},
8592
{
8693
"pointer": "/num_threads",
8794
"type": "int",

‎components/tet_remeshing/wmtk/components/tet_remeshing/tet_remeshing_spec.json‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
"ijk_to_ras",
99
"skip_simplify",
1010
"use_sample_envelope",
11+
"use_tetgen",
1112
"num_threads",
1213
"max_iterations",
1314
"eps_rel",
@@ -77,6 +78,12 @@
7778
"default": false,
7879
"doc": "Use sample envelope instead of exact one."
7980
},
81+
{
82+
"pointer": "/use_tetgen",
83+
"type": "bool",
84+
"default": false,
85+
"doc": "Use tetgen for embedding an image. Potentially faster but could fail."
86+
},
8087
{
8188
"pointer": "/num_threads",
8289
"type": "int",

0 commit comments

Comments
 (0)