integration-of-complex-surfaces
/
ImplicitSurfaceNetwork
mirror of http://118.195.195.192:3000/ZCWang/ImplicitSurfaceNetwork.git
2
0
Fork 1
Code Issues Projects Releases Wiki Activity
Browse Source

fix runtime error; 

improve the code structure of extract_mesh and pair_faces
V2-integral-fix
Zhicheng Wang 4 weeks ago
parent
e0d2500971
commit
024386ea98
6 changed files with 569 additions and 520 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 17
      network_process/interface/fwd_types.hpp
  2. 541
      network_process/src/connect_by_topo/pair_faces.cpp
  3. 2
      network_process/src/connect_by_topo/patch_connectivity.cpp
  4. 197
      network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp
  5. 322
      network_process/src/prim_gen/extract_mesh.cpp
  6. 10
      network_process/src/process.cpp

17
network_process/interface/fwd_types.hpp
View File

@ -29,6 +29,23 @@ static inline uint32_t hash_vert_pos(size_t x, size_t y, size_t z, const scene_b
return static_cast<uint32_t>(z + info.vert_resolution.z() * (y + info.vert_resolution.y() * x));
}
static inline uint32_t pos_hash_max(const scene_bg_mesh_info_t& info)
{
return hash_vert_pos(info.vert_resolution.x() - 1, info.vert_resolution.y() - 1, info.vert_resolution.z() - 1, info);
}
static inline uint32_t hash_increment_x(uint32_t hashed_value, const scene_bg_mesh_info_t& info)
{
return hashed_value + static_cast<uint32_t>(info.vert_resolution.y() * info.vert_resolution.z());
}
static inline uint32_t hash_increment_y(uint32_t hashed_value, const scene_bg_mesh_info_t& info)
{
return hashed_value + static_cast<uint32_t>(info.vert_resolution.z());
}
static inline uint32_t hash_increment_z(uint32_t hashed_value, const scene_bg_mesh_info_t& info) { return hashed_value + 1; }
static inline Eigen::Vector3d get_grid_pos(size_t x, size_t y, size_t z, const scene_bg_mesh_info_t& info)
{
return info.grid_offset + Eigen::Vector3d(x, y, z).cwiseProduct(info.grid_size);

541
network_process/src/connect_by_topo/pair_faces.cpp
View File

@ -1,5 +1,12 @@
#include <functional>
#include <numeric>
#include <connect_by_topo.hpp>
/// TODO: for applying parallelism, take every branch separately:
/// 1. pair in one tet;
/// 2. calculate if on_tet_edge;
/// 2.1. not on_tet_edge, and containing_simplex.size() == 3;
/// 2.2. containing_simplex.size() < 3
void compute_patch_order(const stl_vector_mp<std::array<uint32_t, 4>> &tetrahedrons,
const chain_header_t &character_edge,
const stl_vector_mp<iso_vertex_t> &iso_verts,
@ -28,42 +35,53 @@ void compute_patch_order(const stl_vector_mp<std::array<uint32_t, 4>>
patch_of_face_mapping,
half_patch_adj_list);
} else {
const auto v1 = character_edge.v1;
const auto v2 = character_edge.v2;
const auto &iso_v1 = iso_verts[v1];
const auto &iso_v2 = iso_verts[v2];
bool on_tet_edge = false;
uint32_t vId1, vId2, vId3;
if (iso_v1.header.minimal_simplex_flag == 1 && iso_v2.header.minimal_simplex_flag == 1) {
on_tet_edge = true;
vId1 = iso_v1.simplex_vertex_indices[0];
vId2 = iso_v2.simplex_vertex_indices[0];
} else if (iso_v1.header.minimal_simplex_flag == 2 && iso_v2.header.minimal_simplex_flag == 1) {
vId1 = iso_v1.simplex_vertex_indices[0];
vId2 = iso_v1.simplex_vertex_indices[1];
vId3 = iso_v2.simplex_vertex_indices[0];
on_tet_edge = (vId3 == vId1 || vId3 == vId2);
} else if (iso_v1.header.minimal_simplex_flag == 1 && iso_v2.header.minimal_simplex_flag == 2) {
vId1 = iso_v2.simplex_vertex_indices[0];
vId2 = iso_v2.simplex_vertex_indices[1];
vId3 = iso_v1.simplex_vertex_indices[0];
on_tet_edge = (vId3 == vId1 || vId3 == vId2);
} else if (iso_v1.header.minimal_simplex_flag == 2 && iso_v2.header.minimal_simplex_flag == 2) {
vId1 = iso_v1.simplex_vertex_indices[0];
vId2 = iso_v1.simplex_vertex_indices[1];
bool on_tet_edge{false};
uint32_t vId1, vId2;
{
const auto v1 = character_edge.v1;
const auto v2 = character_edge.v2;
auto iso_v1_ptr = &iso_verts[v1];
auto iso_v2_ptr = &iso_verts[v2];
auto minimal_simplex_count1 = iso_v1_ptr->header.minimal_simplex_flag;
auto minimal_simplex_count2 = iso_v2_ptr->header.minimal_simplex_flag;
if (minimal_simplex_count1 > minimal_simplex_count2) {
std::swap(iso_v1_ptr, iso_v2_ptr);
std::swap(minimal_simplex_count1, minimal_simplex_count2);
}
const auto &iso_v1 = *iso_v1_ptr;
const auto &iso_v2 = *iso_v2_ptr;
const auto &simplex_vertex_indices1 = iso_v1.simplex_vertex_indices;
const auto &simplex_vertex_indices2 = iso_v2.simplex_vertex_indices;
// assume: indices in Iso_Vert::simplex_vertex_indices are sorted
on_tet_edge = (vId1 == iso_v2.simplex_vertex_indices[0] && vId2 == iso_v2.simplex_vertex_indices[1]);
if (minimal_simplex_count1 == 1 && minimal_simplex_count2 == 1) {
on_tet_edge = true;
vId1 = simplex_vertex_indices1.front();
vId2 = simplex_vertex_indices2.front();
} else if (minimal_simplex_count1 == 1 && minimal_simplex_count2 == 2) {
vId1 = simplex_vertex_indices2[0];
vId2 = simplex_vertex_indices2[1];
auto start_iter = simplex_vertex_indices2.begin();
auto end_iter = simplex_vertex_indices2.begin() + minimal_simplex_count2;
auto iter = std::find(start_iter, end_iter, simplex_vertex_indices1.front());
on_tet_edge = iter != end_iter;
} else if (minimal_simplex_count1 == 2 && minimal_simplex_count2 == 2) {
vId1 = simplex_vertex_indices1[0];
vId2 = simplex_vertex_indices1[1];
on_tet_edge = simplex_vertex_indices1 == simplex_vertex_indices2;
}
}
if (on_tet_edge) {
// character_edge lies on tet edge (vId1, vId2)
// tet vertices vId1, vId2 must be degenerate vertices
// find all tets incident to edge (vId1, vId2)
const auto &tet_ids = zero_vertex_to_incident_tet_mapping.at(vId1);
flat_hash_set_mp<uint32_t> incident_tets1(tet_ids.begin(), tet_ids.end());
stl_vector_mp<uint32_t> common_tetIds{};
for (const auto tId : zero_vertex_to_incident_tet_mapping.at(vId2)) {
if (incident_tets1.find(tId) != incident_tets1.end()) { common_tetIds.emplace_back(tId); }
}
const auto &tet_ids1 = zero_vertex_to_incident_tet_mapping.at(vId1);
const auto &tet_ids2 = zero_vertex_to_incident_tet_mapping.at(vId2);
stl_vector_mp<uint32_t> common_tetIds{};
std::set_intersection(tet_ids1.begin(),
tet_ids1.end(),
tet_ids2.begin(),
tet_ids2.end(),
std::back_inserter(common_tetIds));
// pair half-faces
pair_patches_in_tets({vId1, vId2},
common_tetIds,
@ -78,14 +96,17 @@ void compute_patch_order(const stl_vector_mp<std::array<uint32_t, 4>>
} else {
// character_edge lies on a tet boundary face
// assert: containing_tets.size() == 2
const auto tet_id1 = character_edge.volume_indices.front();
const auto tet_id2 = character_edge.volume_indices.back();
flat_hash_set_mp<uint32_t> tet1_vIds(tetrahedrons[tet_id1].begin(), tetrahedrons[tet_id1].end());
stl_vector_mp<uint32_t> common_vIds{};
const auto tet_id1 = character_edge.volume_indices.front();
const auto tet_id2 = character_edge.volume_indices.back();
const auto &vert_ids1 = tetrahedrons[tet_id1];
const auto &vert_ids2 = tetrahedrons[tet_id2];
stl_vector_mp<uint32_t> common_vIds{};
common_vIds.reserve(3);
for (const auto vId : tetrahedrons[tet_id2]) {
if (tet1_vIds.find(vId) != tet1_vIds.end()) common_vIds.emplace_back(vId);
}
std::set_intersection(vert_ids1.begin(),
vert_ids1.end(),
vert_ids2.begin(),
vert_ids2.end(),
std::back_inserter(common_vIds));
// pair half-faces
pair_patches_in_tets(common_vIds,
{tet_id1, tet_id2},
@ -102,6 +123,12 @@ void compute_patch_order(const stl_vector_mp<std::array<uint32_t, 4>>
}
// ===============================================================================================
struct travel_info_t {
// uint32_t iso_face_id{invalid_index};
uint32_t face_id{invalid_index};
int8_t face_sign{0};
};
// compute neighboring pair of half-patches around an iso-edge in a tetrahedron
// half-patch adjacency list : (patch i, 1) <--> 2i, (patch i, -1) <--> 2i+1
void pair_patches_in_one_tet(const arrangement_t &tet_cut_result,
@ -111,85 +138,74 @@ void pair_patches_in_one_tet(const arrangement_t &tet_cut_res
stl_vector_mp<stl_vector_mp<uint32_t>> &half_patch_adj_list)
{
// find tet faces that are incident to the character_edge
stl_vector_mp<bool> is_incident_faces(tet_cut_result.faces.size(), false);
dynamic_bitset_mp<> is_incident_faces(tet_cut_result.faces.size());
// map: tet face id --> iso face id
stl_vector_mp<uint32_t> iso_face_Id_of_face(tet_cut_result.faces.size(), invalid_index);
stl_vector_mp<uint32_t> iso_face_id_of_face(tet_cut_result.faces.size(), invalid_index);
for (const auto &iso_face_id : character_edge.face_indices) {
const auto face_id = iso_faces[iso_face_id].headers[0].local_face_index;
const auto face_id = iso_faces[iso_face_id].headers.front().local_face_index;
is_incident_faces[face_id] = true;
iso_face_Id_of_face[face_id] = iso_face_id;
iso_face_id_of_face[face_id] = iso_face_id;
}
// travel around the edge
stl_vector_mp<bool> visited_cell(tet_cut_result.cells.size(), false);
struct travel_info_t {
uint32_t iso_face_id{invalid_index};
uint32_t face_id{invalid_index};
int8_t face_sign{0};
void clear()
{
iso_face_id = invalid_index;
face_id = invalid_index;
face_sign = 0;
}
void move_update(travel_info_t &&other)
{
*this = std::move(other);
this->face_sign = -this->face_sign;
other.clear();
}
auto get_half_patch_index = [&](const travel_info_t &info) {
const auto iso_face_index = iso_face_id_of_face[info.face_id];
const auto patch_index = patch_of_face_mapping[iso_face_index];
return info.face_sign > 0 ? 2 * patch_index : 2 * patch_index + 1;
};
auto travel_func = [&](uint32_t start_face_index, int8_t start_sign) {
const auto &faces = tet_cut_result.faces;
travel_info_t info{iso_faces[start_face_index].headers.front().local_face_index, start_sign};
travel_info_t other_info{};
auto cell_id = start_sign > 0 ? faces[info.face_id].positive_cell : faces[info.face_id].negative_cell;
auto travel_func = [&](uint32_t &cell_id, travel_info_t &info1, travel_info_t &info2) {
while (cell_id != invalid_index && !visited_cell[cell_id]) {
visited_cell[cell_id] = true;
const auto &cell_faces = tet_cut_result.cells[cell_id].faces;
visited_cell[cell_id] = true;
// find next face
for (const auto &fId : tet_cut_result.cells[cell_id].faces) {
if (is_incident_faces[fId] && fId != info1.face_id) { info2.face_id = fId; }
}
if (info2.face_id == invalid_index) {
auto iter = std::find_if(cell_faces.begin(), cell_faces.end(), [&](uint32_t fId) {
return is_incident_faces[fId] && fId != info.face_id;
});
if (iter == cell_faces.end()) {
// next face is not found
break;
} else {
// get sign of face2 and find next cell
if (tet_cut_result.faces[info2.face_id].positive_cell == cell_id) {
cell_id = tet_cut_result.faces[info2.face_id].negative_cell;
info2.face_sign = 1;
if (faces[*iter].positive_cell == cell_id) {
cell_id = faces[*iter].negative_cell;
other_info.face_sign = 1;
} else {
cell_id = tet_cut_result.faces[info2.face_id].positive_cell;
info2.face_sign = -1;
cell_id = faces[*iter].positive_cell;
other_info.face_sign = -1;
}
other_info.face_id = *iter;
// add (face1, face2) to the list of face pairs
info2.iso_face_id = iso_face_Id_of_face[info2.face_id];
const auto half_patch_index1 = info1.face_sign > 0 ? 2 * patch_of_face_mapping[info1.iso_face_id]
: 2 * patch_of_face_mapping[info1.iso_face_id] + 1;
const auto half_patch_index2 = info2.face_sign > 0 ? 2 * patch_of_face_mapping[info2.iso_face_id]
: 2 * patch_of_face_mapping[info2.iso_face_id] + 1;
const auto half_patch_index1 = get_half_patch_index(info);
const auto half_patch_index2 = get_half_patch_index(other_info);
half_patch_adj_list[half_patch_index1].emplace_back(half_patch_index2);
half_patch_adj_list[half_patch_index2].emplace_back(half_patch_index1);
// update face1 and clear face2
info1.move_update(std::move(info2));
info = other_info;
info.face_sign = -info.face_sign;
}
}
};
const auto character_face_index = character_edge.face_indices.front();
const auto character_face_index = character_edge.face_indices.front();
// start from the first incident face, and travel to its positive cell
travel_info_t info1{character_face_index, iso_faces[character_face_index].headers[0].local_face_index, 1};
travel_info_t info2{};
auto cell_id = tet_cut_result.faces[info1.face_id].positive_cell;
travel_func(cell_id, info1, info2);
travel_func(character_face_index, 1);
// travel in a different direction
info1 = {character_face_index, iso_faces[character_face_index].headers[0].local_face_index, -1};
info2.clear();
cell_id = tet_cut_result.faces[info1.face_id].negative_cell;
travel_func(cell_id, info1, info2);
travel_func(character_face_index, -1);
}
// ===============================================================================================
struct tet_face_travel_info : face_header_t {
int8_t orientation{};
};
const std::array tet_boundary_faces = {0u, 1u, 2u, 3u};
// compute neighboring pair of half-patches around an iso-edge in multiple tetrahedrons
// half-patch adjacency list : (patch i, 1) <--> 2i, (patch i, -1) <--> 2i+1
void pair_patches_in_tets(const stl_vector_mp<uint32_t> &containing_simplex,
@ -206,38 +222,30 @@ void pair_patches_in_tets(const stl_vector_mp<uint32_t> &containi
//// pre-processing
// collect all iso-faces incident to the iso-edge
// map: (tet_id, tet_face_id) -> iso_face_id
flat_hash_map_mp<face_header_t, uint32_t> iso_face_Id_of_face{};
flat_hash_map_mp<face_header_t, uint32_t> iso_face_id_of_face{};
for (const auto &iso_face_id : character_edge.face_indices) {
for (const auto &t : iso_faces[iso_face_id].headers) { iso_face_Id_of_face[t] = iso_face_id; }
for (const auto &t : iso_faces[iso_face_id].headers)
iso_face_id_of_face.lazy_emplace(t, [&](const auto &ctor) { ctor(t, iso_face_id); });
}
// find identical tet boundary planes incident to iso-edge
// map: (tet_Id, tet_plane_Id) -> (oppo_tet_Id, oppo_tet_plane_Id)
flat_hash_map_mp<face_header_t, face_header_t> identical_tet_planes{};
if (containing_simplex.size() == 3) {
// character_edge contained in a tet boundary triangle
const auto vId1 = containing_simplex[0];
const auto vId2 = containing_simplex[1];
const auto vId3 = containing_simplex[2];
const auto tId1 = containing_tetIds[0];
const auto tId2 = containing_tetIds[1];
const auto &tet1 = tetrahedrons[tId1];
const auto &tet2 = tetrahedrons[tId2];
uint32_t pId1, pId2;
uint32_t vId;
for (uint32_t i = 0; i < 4; ++i) {
vId = tet1[i];
if (vId != vId1 && vId != vId2 && vId != vId3) {
pId1 = i;
break;
}
}
for (uint32_t i = 0; i < 4; ++i) {
vId = tet2[i];
if (vId != vId1 && vId != vId2 && vId != vId3) {
pId2 = i;
break;
}
}
const auto tId1 = containing_tetIds[0];
const auto tId2 = containing_tetIds[1];
const auto &tet1 = tetrahedrons[tId1];
const auto &tet2 = tetrahedrons[tId2];
auto iter1 = std::find_if(tet1.begin(), tet1.end(), [&](uint32_t vertex_index) {
return std::find(containing_simplex.begin(), containing_simplex.end(), vertex_index) == containing_simplex.end();
});
assert(iter1 != tet1.end());
uint32_t pId1 = std::distance(tet1.begin(), iter1);
auto iter2 = std::find_if(tet2.begin(), tet2.end(), [&](uint32_t vertex_index) {
return std::find(containing_simplex.begin(), containing_simplex.end(), vertex_index) == containing_simplex.end();
});
assert(iter2 != tet2.end());
uint32_t pId2 = std::distance(tet2.begin(), iter2);
// (tId1, pId1) and (tId2, pId2)
identical_tet_planes[{tId1, pId1}] = {tId2, pId2};
identical_tet_planes[{tId2, pId2}] = {tId1, pId1};
@ -266,217 +274,216 @@ void pair_patches_in_tets(const stl_vector_mp<uint32_t> &containi
// find identical faces (tet_id, tet_face_id) on tet boundary planes incident to iso-edge
// map: (tet_Id, tet_face_Id) -> (oppo_tet_Id, oppo_tet_face_Id)
flat_hash_map_mp<face_header_t, face_header_t> opposite_face{};
// map: (tet_Id, tet_vert_Id) -> boundary vert Id
flat_hash_map_mp<face_header_t, uint32_t> boundary_vert_id{};
uint32_t num_boundary_vert = 0;
// hash table for vertices on the boundary of tetrahedron
flat_hash_map_mp<uint32_t, uint32_t> vert_on_tetVert{};
flat_hash_map_mp<pod_key_t<3>, uint32_t> vert_on_tetEdge{};
flat_hash_map_mp<pod_key_t<5>, uint32_t> vert_on_tetFace{};
// keys of vertices on the boundary of tetrahedron
stl_vector_mp<uint32_t> tet_vert_keys{};
stl_vector_mp<pod_key_t<3>> tet_edge_keys{};
stl_vector_mp<pod_key_t<5>> tet_face_keys{};
// hash table for faces on the boundary of tetrahedron
// map: (i,j,k) -> (tet_Id, tet_face_Id)
flat_hash_map_mp<pod_key_t<3>, face_header_t> face_on_tetFace{};
// auxiliary data
stl_vector_mp<bool> is_boundary_vert{};
stl_vector_mp<bool> is_boundary_face{};
flat_hash_set_mp<uint32_t> boundary_vertices{};
flat_hash_set_mp<uint32_t> boundary_faces{};
// map: local vert index --> boundary vert index
stl_vector_mp<uint32_t> boundary_vId_of_vert{};
stl_vector_mp<uint32_t> face_verts{};
pod_key_t<3> key3;
pod_key_t<5> key5;
std::array<bool, 4> used_pId;
std::array<uint32_t, 2> vIds2;
std::array<uint32_t, 3> vIds3;
std::array<uint32_t, 3> implicit_pIds;
std::array<uint32_t, 3> boundary_pIds;
for (const auto i : containing_tetIds) {
std::array<uint32_t, 3> not_boundary_pIds;
for (const auto tet_index : containing_tetIds) {
// non-empty tet i
const auto &arrangement = tetrahedron_arrangements[i];
const auto &arrangement = tetrahedron_arrangements[tet_index];
const auto &vertices = arrangement.vertices;
const auto &faces = arrangement.faces;
auto start_index = tetrahedron_active_subface_start_index[i];
auto num_func = tetrahedron_active_subface_start_index[i + 1] - start_index;
auto start_index = tetrahedron_active_subface_start_index[tet_index];
auto num_func = tetrahedron_active_subface_start_index[tet_index + 1] - start_index;
// find vertices and faces on tet boundary incident to iso-edge
is_boundary_vert.assign(arrangement.vertices.size(), false);
is_boundary_face.clear();
is_boundary_face.reserve(faces.size());
for (const auto &face : faces) {
is_boundary_face.emplace_back(false);
boundary_vertices.clear();
boundary_faces.clear();
for (size_t i = 0; i < faces.size(); ++i) {
const auto &face = faces[i];
if (face.supporting_plane < 4
&& identical_tet_planes.find({i, face.supporting_plane}) != identical_tet_planes.end()) {
is_boundary_face.back() = true;
for (const auto &vid : face.vertices) { is_boundary_vert[vid] = true; }
&& identical_tet_planes.find({tet_index, face.supporting_plane}) != identical_tet_planes.end()) {
boundary_faces.emplace(i);
for (const auto &vertex_index : face.vertices) boundary_vertices.emplace(vertex_index);
}
}
// map: local vert index --> boundary vert index
boundary_vId_of_vert.clear();
boundary_vId_of_vert.reserve(vertices.size());
// create boundary vertices
for (uint32_t j = 0; j < vertices.size(); j++) {
boundary_vId_of_vert.emplace_back(invalid_index);
if (is_boundary_vert[j]) {
uint32_t num_boundary_planes = 0;
uint32_t num_impl_planes = 0;
const auto &vertex = vertices[j];
// vertex.size() == 3
for (const auto &vert_plane : vertex) {
if (vert_plane >= 4) { // plane 0,1,2,3 are tet boundaries
implicit_pIds[num_impl_planes++] = active_subface_indices[vert_plane - 4 + start_index];
} else {
boundary_pIds[num_boundary_planes++] = vert_plane;
}
for (const auto &vertex_index : boundary_vertices) {
const auto &vertex = vertices[vertex_index];
auto &local_vert_to_boundary_vert = boundary_vId_of_vert.emplace_back(invalid_index);
uint32_t num_boundary_planes{};
uint32_t num_impl_planes{};
// vertex.size() == 3
for (const auto &vert_plane : vertex) {
if (vert_plane >= 4) { // plane 0,1,2,3 are tet boundaries
implicit_pIds[num_impl_planes++] = active_subface_indices[vert_plane - 4 + start_index];
} else {
boundary_pIds[num_boundary_planes++] = vert_plane;
}
switch (num_boundary_planes) {
case 2: // on tet edge
{
used_pId.fill(false);
used_pId[boundary_pIds[0]] = true;
used_pId[boundary_pIds[1]] = true;
uint32_t num_vIds = 0;
for (uint32_t k = 0; k < 4; k++) {
if (!used_pId[k]) vIds2[num_vIds++] = tetrahedrons[i][k];
}
uint32_t vId1 = vIds2[0];
uint32_t vId2 = vIds2[1];
if (vId1 > vId2) std::swap(vId1, vId2);
key3[0] = vId1;
key3[1] = vId2;
key3[2] = implicit_pIds[0];
auto iter_inserted = vert_on_tetEdge.try_emplace(key3, num_boundary_vert);
if (iter_inserted.second) { num_boundary_vert++; }
boundary_vId_of_vert.back() = iter_inserted.first->second;
break;
}
case 1: // on tet face
{
uint32_t pId = boundary_pIds[0];
uint32_t num_vIds = 0;
for (uint32_t k = 0; k < 4; k++) {
if (k != pId) vIds3[num_vIds++] = tetrahedrons[i][k];
}
std::sort(vIds3.begin(), vIds3.end());
key5 = {vIds3[0], vIds3[1], vIds3[2], implicit_pIds[0], implicit_pIds[1]};
auto iter_inserted = vert_on_tetFace.try_emplace(key5, num_boundary_vert);
if (iter_inserted.second) { num_boundary_vert++; }
boundary_vId_of_vert.back() = iter_inserted.first->second;
break;
}
case 3: // on tet vertex
{
used_pId.fill(false);
for (const auto &pId : boundary_pIds) { used_pId[pId] = true; }
uint32_t vId;
for (uint32_t k = 0; k < 4; k++) {
if (!used_pId[k]) {
vId = k;
break;
}
}
auto key = tetrahedrons[i][vId];
auto iter_inserted = vert_on_tetVert.try_emplace(key, num_boundary_vert);
if (iter_inserted.second) { num_boundary_vert++; }
boundary_vId_of_vert.back() = iter_inserted.first->second;
break;
}
default: break;
}
std::sort(boundary_pIds.begin(), boundary_pIds.begin() + num_boundary_planes);
std::set_difference(tet_boundary_faces.begin(),
tet_boundary_faces.end(),
boundary_pIds.begin(),
boundary_pIds.begin() + num_impl_planes,
not_boundary_pIds.begin());
const auto &tet_vertices = tetrahedrons[tet_index];
if (num_boundary_planes == 1) { // on tet face
key5 = {tet_vertices[not_boundary_pIds[0]],
tet_vertices[not_boundary_pIds[1]],
tet_vertices[not_boundary_pIds[2]],
implicit_pIds[0],
implicit_pIds[1]};
auto iter = std::find(tet_face_keys.begin(), tet_face_keys.end(), key5);
if (iter != tet_face_keys.end())
local_vert_to_boundary_vert = std::distance(tet_face_keys.begin(), iter);
else {
local_vert_to_boundary_vert = tet_face_keys.size();
tet_face_keys.emplace_back(key5);
}
} else if (num_boundary_planes == 2) { // on tet edge
key3 = {tet_vertices[not_boundary_pIds[0]], tet_vertices[not_boundary_pIds[1]], implicit_pIds[0]};
auto iter = std::find(tet_edge_keys.begin(), tet_edge_keys.end(), key3);
if (iter != tet_edge_keys.end())
local_vert_to_boundary_vert = std::distance(tet_edge_keys.begin(), iter);
else {
local_vert_to_boundary_vert = tet_edge_keys.size();
tet_edge_keys.emplace_back(key3);
}
} else if (num_boundary_planes == 3) { // on tet vertex
auto key = tet_vertices[not_boundary_pIds[0]];
auto iter = std::find(tet_vert_keys.begin(), tet_vert_keys.end(), key);
if (iter != tet_vert_keys.end())
local_vert_to_boundary_vert = std::distance(tet_vert_keys.begin(), iter);
else {
local_vert_to_boundary_vert = tet_vert_keys.size();
tet_vert_keys.emplace_back(key);
}
}
}
// pair boundary faces
for (uint32_t j = 0; j < faces.size(); ++j) {
if (is_boundary_face[j]) {
face_verts.clear();
for (auto vId : faces[j].vertices) face_verts.emplace_back(boundary_vId_of_vert[vId]);
compute_iso_face_key(face_verts, key3);
auto iter_inserted = face_on_tetFace.try_emplace(key3, face_header_t{i, j});
if (!iter_inserted.second) {
// face inserted before, pair the two faces
opposite_face[iter_inserted.first->second] = {i, j};
opposite_face[{i, j}] = iter_inserted.first->second;
}
for (const auto &face_index : boundary_faces) {
const auto &face = faces[face_index];
face_verts.clear();
for (auto vId : face.vertices) face_verts.emplace_back(boundary_vId_of_vert[vId]);
compute_iso_face_key(face_verts, key3);
face_header_t face_header{tet_index, face_index};
auto iter_inserted = face_on_tetFace.try_emplace(key3, face_header);
if (!iter_inserted.second) {
// face inserted before, pair the two faces
const auto &other_face_header = iter_inserted.first->second;
opposite_face[face_header] = other_face_header;
opposite_face[other_face_header] = face_header;
}
}
}
// find the half-iso-face of the local face in tet with given orientation
// the orientation of an iso-face is defined by the smallest-index implicit function passing the iso-face
auto get_half_iso_face = [&](face_header_t tet_face, int8_t orient, uint32_t &iso_face_id, int8_t &iso_orient) {
iso_face_id = iso_face_Id_of_face[tet_face];
const auto &cell_complex = tetrahedron_arrangements[tet_face.volume_index];
const auto &faces = cell_complex.faces;
auto supp_pId = faces[tet_face.local_face_index].supporting_plane;
if (supp_pId > 3) { // plane 0,1,2,3 are tet boundary planes
// supporting plane is not a tet boundary plane
// assume: supporting plane is the iso-surface of smallest-index implicit function
iso_orient = orient;
} else {
auto get_half_iso_face = [&](tet_face_travel_info tet_info) {
const auto iso_face_index = iso_face_id_of_face[tet_info];
const auto &cell_complex = tetrahedron_arrangements[tet_info.volume_index];
const auto &faces = cell_complex.faces;
auto supp_pId = faces[tet_info.local_face_index].supporting_plane;
int8_t iso_face_orient = tet_info.orientation;
// if (supp_pId > 3) { // plane 0,1,2,3 are tet boundary planes
// // supporting plane is not a tet boundary plane
// // assume: supporting plane is the iso-surface of smallest-index implicit function
// // note: this is because when the second or more implicit function's iso-surface is constructing, it's an split
// of
// // the first iso-surface
// iso_face_orient = tet_info.orientation;
// } else {
// // supporting plane is a tet boundary plane, must be duplicate planes
// const auto uId = cell_complex.unique_plane_indices[supp_pId];
// // find the smallest-index non-boundary plane
// uint32_t min_pId{invalid_index};
// for (auto pId : cell_complex.unique_planes[uId]) {
// if (pId > 3 && pId < min_pId) { min_pId = pId; }
// }
// // orient is the orientation of supporting plane
// // flip orientation if smallest-index non-boundary plane has different orientation
// if (cell_complex.unique_plane_orientations[min_pId] != cell_complex.unique_plane_orientations[supp_pId]) {
// iso_face_orient = -tet_info.orientation;
// }
// }
if (supp_pId <= 3) {
// supporting plane is a tet boundary plane, must be duplicate planes
const auto uid = cell_complex.unique_plane_indices[supp_pId];
const auto uId = cell_complex.unique_plane_indices[supp_pId];
// find the smallest-index non-boundary plane
uint32_t min_pId = std::numeric_limits<uint32_t>::max();
for (auto pId : cell_complex.unique_planes[uid]) {
if (pId > 3 && pId < min_pId) { min_pId = pId; }
}
uint32_t min_pId = std::transform_reduce(
cell_complex.unique_planes[uId].begin(),
cell_complex.unique_planes[uId].end(),
invalid_index,
[](uint32_t lhs, uint32_t rhs) { return std::min(lhs, rhs); },
[](uint32_t pId) { return pId > 3 ? pId : invalid_index; });
// for (auto pId : cell_complex.unique_planes[uId]) {
// if (pId > 3 && pId < min_pId) { min_pId = pId; }
// }
// orient is the orientation of supporting plane
// flip orientation if smallest-index non-boundary plane has different orientation
if (cell_complex.unique_plane_orientations[min_pId] != cell_complex.unique_plane_orientations[supp_pId]) {
iso_orient = -orient;
} else {
iso_orient = orient;
iso_face_orient = -tet_info.orientation;
}
}
return iso_face_orient > 0 ? 2 * patch_of_face_mapping[iso_face_index] : 2 * patch_of_face_mapping[iso_face_index] + 1;
};
// pair (tet_id, tet_face_id)
auto find_next =
[&](face_header_t face, int8_t orient, face_header_t &face_next, int8_t &orient_next, auto &&find_next) -> void {
const auto tet_id = face.volume_index;
const auto tet_face_id = face.local_face_index;
auto find_next = [&](const tet_face_travel_info &info_cur, tet_face_travel_info &info_next, auto &&find_next) -> void {
const auto tet_id = info_cur.volume_index;
const auto tet_face_id = info_cur.local_face_index;
// non-empty tet
const auto &cell_complex = tetrahedron_arrangements[tet_id];
uint32_t cell_id =
(orient == 1 ? cell_complex.faces[tet_face_id].positive_cell : cell_complex.faces[tet_face_id].negative_cell);
const auto &complex = tetrahedron_arrangements[tet_id];
const auto &faces = complex.faces;
uint32_t cell_id = (info_cur.orientation == 1 ? faces[tet_face_id].positive_cell : faces[tet_face_id].negative_cell);
if (cell_id != invalid_index) {
for (auto fi : cell_complex.cells[cell_id].faces) {
if (fi != tet_face_id
&& (iso_face_Id_of_face.find({tet_id, fi}) != iso_face_Id_of_face.end()
|| opposite_face.find({tet_id, fi}) != opposite_face.end())) {
face_next = {tet_id, fi};
break;
}
}
orient_next = cell_complex.faces[face_next.local_face_index].positive_cell == cell_id ? 1 : -1;
const auto &cell_faces = complex.cells[cell_id].faces;
auto iter = std::find_if(cell_faces.begin(), cell_faces.end(), [&](uint32_t face_index) {
const face_header_t other_header{tet_id, face_index};
return face_index != tet_face_id && iso_face_id_of_face.find(other_header) != iso_face_id_of_face.end()
|| opposite_face.find(other_header) != opposite_face.end();
});
assert(iter != cell_faces.end());
info_next.volume_index = tet_id;
info_next.local_face_index = *iter;
info_next.orientation = faces[info_next.local_face_index].positive_cell == cell_id ? 1 : -1;
} else {
// cell is None, so the face lies on tet boundary
find_next(opposite_face[face], 1, face_next, orient_next, find_next);
tet_face_travel_info other_info_cur{opposite_face[info_cur], 1};
find_next(other_info_cur, info_next, find_next);
}
};
//// face pairing algorithm
const auto iso_face_0 = character_edge.face_indices.front();
const auto iso_face_0 = character_edge.face_indices.front();
// pair (tet_id, tet_face_id)
auto face_curr = iso_faces[iso_face_0].headers[0];
int8_t orient_curr = 1; // orientation: 1 for Positive, -1 for Negative
face_header_t face_next;
int8_t orient_next, iso_orient_curr, iso_orient_next;
uint32_t iso_face_curr = iso_face_0;
uint32_t iso_face_next = invalid_index;
while (iso_face_next != iso_face_0) {
find_next(face_curr, orient_curr, face_next, orient_next, find_next);
while (iso_face_Id_of_face.find(face_next) == iso_face_Id_of_face.end()) {
find_next(face_next, -orient_next, face_next, orient_next, find_next);
tet_face_travel_info tet_info_cur{iso_faces[iso_face_0].headers.front(), 1}, tet_info_next;
uint32_t iso_face_index_next{invalid_index};
while (iso_face_index_next != iso_face_0) {
find_next(tet_info_cur, tet_info_next, find_next);
while (iso_face_id_of_face.find(tet_info_next) == iso_face_id_of_face.end()) {
tet_info_next.orientation = -tet_info_next.orientation;
find_next(tet_info_next, tet_info_next, find_next);
}
get_half_iso_face(face_curr, orient_curr, iso_face_curr, iso_orient_curr);
get_half_iso_face(face_next, orient_next, iso_face_next, iso_orient_next);
//
const auto half_patch_index1 =
iso_orient_curr > 0 ? 2 * patch_of_face_mapping[iso_face_curr] : 2 * patch_of_face_mapping[iso_face_curr] + 1;
const auto half_patch_index2 =
iso_orient_next > 0 ? 2 * patch_of_face_mapping[iso_face_next] : 2 * patch_of_face_mapping[iso_face_next] + 1;
const auto half_patch_index1 = get_half_iso_face(tet_info_cur);
const auto half_patch_index2 = get_half_iso_face(tet_info_next);
half_patch_adj_list[half_patch_index1].emplace_back(half_patch_index2);
half_patch_adj_list[half_patch_index2].emplace_back(half_patch_index1);
//
face_curr = face_next;
orient_curr = -orient_next;
iso_face_index_next = iso_face_id_of_face[tet_info_cur];
tet_info_cur = tet_info_next;
tet_info_cur.orientation = -tet_info_cur.orientation;
}
}

2
network_process/src/connect_by_topo/patch_connectivity.cpp
View File

@ -13,7 +13,7 @@ void compute_patch_edges(const stl_vector_mp<polygon_face_t>& patch_face
edges_of_face.reserve(patch_faces.size());
uint32_t num_iso_edge{};
// map: (v1, v2) -> iso-edge index
flat_hash_map_mp<std::pair<uint32_t, uint32_t>, uint32_t> edge_id{};
// flat_hash_map_mp<std::pair<uint32_t, uint32_t>, uint32_t> edge_id{};
for (uint32_t i = 0; i < patch_faces.size(); i++) {
auto& face = patch_faces[i];
const uint32_t num_edge = static_cast<uint32_t>(face.vertex_indices.size());

197
network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp
View File

@ -21,78 +21,139 @@ void build_tetrahedron_and_adjacency(const scene_bg_mesh_info_t&
flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& reverse_vertex_adjacency,
btree_map_mp<uint32_t, stl_vector_mp<uint32_t>>& vertex_to_tet_mapping)
{
auto max_pos = scene_bg_mesh_info.vert_resolution.template cast<size_t>() - Eigen::Vector<size_t, 3>::Ones();
tetrahedrons.reserve((max_pos.array() - 1).prod() * 5);
scene_bg_mesh_info_t culled_scene_bg_mesh_info = scene_bg_mesh_info;
culled_scene_bg_mesh_info.vert_resolution -= Eigen::Vector<uint64_t, 3>::Ones();
tetrahedrons.reserve((culled_scene_bg_mesh_info.vert_resolution.array() - 1).prod() * 5);
for (size_t x = 0; x < max_pos.x(); ++x) {
for (size_t y = 0; y < max_pos.y(); ++y) {
for (size_t z = 0; z < max_pos.z(); ++z) {
uint32_t v0 = hash_vert_pos(x, y, z, scene_bg_mesh_info);
uint32_t v1 = hash_vert_pos(x + 1, y, z, scene_bg_mesh_info);
uint32_t v2 = hash_vert_pos(x + 1, y + 1, z, scene_bg_mesh_info);
uint32_t v3 = hash_vert_pos(x, y + 1, z, scene_bg_mesh_info);
uint32_t v4 = hash_vert_pos(x, y, z + 1, scene_bg_mesh_info);
uint32_t v5 = hash_vert_pos(x + 1, y, z + 1, scene_bg_mesh_info);
uint32_t v6 = hash_vert_pos(x + 1, y + 1, z + 1, scene_bg_mesh_info);
uint32_t v7 = hash_vert_pos(x, y + 1, z + 1, scene_bg_mesh_info);
auto insert_or_compare_vertex_adjacency = [&](uint32_t src, uint32_t dst) {
reverse_vertex_adjacency[dst].reserve(reverse_vertex_adjacency[dst].size() + 3);
auto [iter, is_new] = vertex_lexigraphical_adjacency.try_emplace(src, dst);
if (!is_new) iter->second = std::min(iter->second, dst);
reverse_vertex_adjacency[dst].emplace_back(src);
};
for (size_t i = 0; i < pos_hash_max(culled_scene_bg_mesh_info); i += 2) {
auto [x, y, z] = get_grid_pos(i, culled_scene_bg_mesh_info);
auto v0 = hash_vert_pos(x, y, z, scene_bg_mesh_info);
auto v1 = hash_increment_x(v0, scene_bg_mesh_info);
auto v2 = hash_increment_y(v1, scene_bg_mesh_info);
auto v3 = hash_increment_y(v0, scene_bg_mesh_info);
auto v4 = hash_increment_z(v0, scene_bg_mesh_info);
auto v5 = hash_increment_x(v4, scene_bg_mesh_info);
auto v6 = hash_increment_y(v5, scene_bg_mesh_info);
auto v7 = hash_increment_y(v4, scene_bg_mesh_info);
tetrahedrons.emplace_back(std::array{v1, v3, v4, v6});
tetrahedrons.emplace_back(std::array{v3, v4, v6, v7});
tetrahedrons.emplace_back(std::array{v0, v1, v3, v4});
tetrahedrons.emplace_back(std::array{v1, v2, v3, v6});
tetrahedrons.emplace_back(std::array{v1, v4, v5, v6});
insert_or_compare_vertex_adjacency(v1, v0);
insert_or_compare_vertex_adjacency(v2, v3);
insert_or_compare_vertex_adjacency(v3, v0);
insert_or_compare_vertex_adjacency(v4, v0);
insert_or_compare_vertex_adjacency(v5, v4);
insert_or_compare_vertex_adjacency(v6, v4);
insert_or_compare_vertex_adjacency(v7, v4);
}
if ((x + y + z) % 2 == 0) {
tetrahedrons.emplace_back(std::array{v4, v6, v1, v3});
tetrahedrons.emplace_back(std::array{v6, v3, v4, v7});
tetrahedrons.emplace_back(std::array{v1, v3, v0, v4});
tetrahedrons.emplace_back(std::array{v3, v1, v2, v6});
tetrahedrons.emplace_back(std::array{v4, v1, v6, v5});
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v6, v1, v3},
v4);
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v6, v3, v7},
v4);
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v1, v3, v4},
v0);
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v1, v2, v6},
v3);
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v1, v6, v5},
v4);
} else {
tetrahedrons.emplace_back(std::array{v7, v0, v2, v5});
tetrahedrons.emplace_back(std::array{v2, v3, v0, v7});
tetrahedrons.emplace_back(std::array{v5, v7, v0, v4});
tetrahedrons.emplace_back(std::array{v7, v2, v6, v5});
tetrahedrons.emplace_back(std::array{v0, v1, v2, v5});
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v7, v2, v5},
v0);
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v2, v3, v7},
v0);
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v5, v7, v4},
v0);
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v2, v6, v5},
v7);
insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
reverse_vertex_adjacency,
{v1, v2, v5},
v0);
}
}
}
for (size_t i = 1; i < pos_hash_max(culled_scene_bg_mesh_info); i += 2) {
auto [x, y, z] = get_grid_pos(i, culled_scene_bg_mesh_info);
auto v0 = hash_vert_pos(x, y, z, scene_bg_mesh_info);
auto v1 = hash_increment_x(v0, scene_bg_mesh_info);
auto v2 = hash_increment_y(v1, scene_bg_mesh_info);
auto v3 = hash_increment_y(v0, scene_bg_mesh_info);
auto v4 = hash_increment_z(v0, scene_bg_mesh_info);
auto v5 = hash_increment_x(v4, scene_bg_mesh_info);
auto v6 = hash_increment_y(v5, scene_bg_mesh_info);
auto v7 = hash_increment_y(v4, scene_bg_mesh_info);
tetrahedrons.emplace_back(std::array{v0, v2, v5, v7});
tetrahedrons.emplace_back(std::array{v0, v2, v3, v7});
tetrahedrons.emplace_back(std::array{v0, v4, v5, v7});
tetrahedrons.emplace_back(std::array{v2, v5, v6, v7});
tetrahedrons.emplace_back(std::array{v0, v1, v2, v5});
insert_or_compare_vertex_adjacency(v1, v0);
insert_or_compare_vertex_adjacency(v2, v0);
insert_or_compare_vertex_adjacency(v3, v0);
insert_or_compare_vertex_adjacency(v4, v0);
insert_or_compare_vertex_adjacency(v5, v0);
insert_or_compare_vertex_adjacency(v6, v7);
insert_or_compare_vertex_adjacency(v7, v0);
}
// for (size_t x = 0; x < max_pos.x(); ++x) {
// for (size_t y = 0; y < max_pos.y(); ++y) {
// for (size_t z = 0; z < max_pos.z(); ++z) {
// uint32_t v0 = hash_vert_pos(x, y, z, scene_bg_mesh_info);
// uint32_t v1 = hash_vert_pos(x + 1, y, z, scene_bg_mesh_info);
// uint32_t v2 = hash_vert_pos(x + 1, y + 1, z, scene_bg_mesh_info);
// uint32_t v3 = hash_vert_pos(x, y + 1, z, scene_bg_mesh_info);
// uint32_t v4 = hash_vert_pos(x, y, z + 1, scene_bg_mesh_info);
// uint32_t v5 = hash_vert_pos(x + 1, y, z + 1, scene_bg_mesh_info);
// uint32_t v6 = hash_vert_pos(x + 1, y + 1, z + 1, scene_bg_mesh_info);
// uint32_t v7 = hash_vert_pos(x, y + 1, z + 1, scene_bg_mesh_info);
// if ((x + y + z) % 2 == 0) {
// tetrahedrons.emplace_back(std::array{v4, v6, v1, v3});
// tetrahedrons.emplace_back(std::array{v6, v3, v4, v7});
// tetrahedrons.emplace_back(std::array{v1, v3, v0, v4});
// tetrahedrons.emplace_back(std::array{v3, v1, v2, v6});
// tetrahedrons.emplace_back(std::array{v4, v1, v6, v5});
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v6, v1, v3},
// v4);
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v6, v3, v7},
// v4);
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v1, v3, v4},
// v0);
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v1, v2, v6},
// v3);
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v1, v6, v5},
// v4);
// } else {
// tetrahedrons.emplace_back(std::array{v7, v0, v2, v5});
// tetrahedrons.emplace_back(std::array{v2, v3, v0, v7});
// tetrahedrons.emplace_back(std::array{v5, v7, v0, v4});
// tetrahedrons.emplace_back(std::array{v7, v2, v6, v5});
// tetrahedrons.emplace_back(std::array{v0, v1, v2, v5});
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v7, v2, v5},
// v0);
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v2, v3, v7},
// v0);
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v5, v7, v4},
// v0);
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v2, v6, v5},
// v7);
// insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
// reverse_vertex_adjacency,
// {v1, v2, v5},
// v0);
// }
// }
// }
// }
for (auto i = 0; i < tetrahedrons.size(); ++i) {
vertex_to_tet_mapping[tetrahedrons[i][0]].emplace_back(i);
vertex_to_tet_mapping[tetrahedrons[i][1]].emplace_back(i);

322
network_process/src/prim_gen/extract_mesh.cpp
View File

@ -1,6 +1,8 @@
#include <prim_gen.hpp>
#include <helper.hpp>
const std::array tet_boundary_faces = {0u, 1u, 2u, 3u};
void extract_iso_mesh(const std::array<uint32_t, 3>& tet_active_subface_counts,
const stl_vector_mp<arrangement_t>& tetrahedron_arrangements,
const scene_bg_mesh_info_t& scene_bg_mesh_info,
@ -32,225 +34,179 @@ void extract_iso_mesh(const std::array<uint32_t, 3>& tet_active_s
flat_hash_map_mp<pod_key_t<3>, uint32_t> face_on_tetFace{};
//
stl_vector_mp<bool> is_iso_vert{};
is_iso_vert.reserve(8);
stl_vector_mp<bool> is_iso_face{};
is_iso_face.reserve(9);
stl_vector_mp<uint32_t> iso_vId_of_vert{};
iso_vId_of_vert.reserve(8);
stl_vector_mp<uint32_t> face_verts{};
flat_hash_set_mp<uint32_t> tet_iso_verts{};
flat_hash_set_mp<uint32_t> tet_iso_faces{};
// map: local vert index --> iso-vert index
stl_vector_mp<uint32_t> iso_vId_of_vert{};
stl_vector_mp<uint32_t> face_verts{};
face_verts.reserve(4);
pod_key_t<3> key3;
pod_key_t<5> key5;
std::array<bool, 4> used_pId;
std::array<uint32_t, 2> vIds2;
std::array<uint32_t, 3> vIds3;
std::array<uint16_t, 3> implicit_pIds;
std::array<uint32_t, 3> boundary_pIds;
std::array<uint32_t, 3> not_boundary_pIds;
//
for (uint32_t i = 0; i < tetrahedrons.size(); i++) {
const auto& arrangement = tetrahedron_arrangements[i];
for (uint32_t tet_index = 0; tet_index < tetrahedrons.size(); tet_index++) {
const auto& arrangement = tetrahedron_arrangements[tet_index];
const auto& vertices = arrangement.vertices;
const auto& faces = arrangement.faces;
auto start_index = tetrahedron_active_subface_start_index[i];
auto start_index = tetrahedron_active_subface_start_index[tet_index];
// find vertices and faces on isosurface
is_iso_vert.assign(vertices.size(), false);
is_iso_face.clear();
is_iso_face.reserve(faces.size());
tet_iso_verts.clear();
tet_iso_verts.reserve(vertices.size());
tet_iso_faces.clear();
tet_iso_faces.reserve(faces.size());
if (arrangement.unique_planes.empty()) { // all planes are unique
for (const auto& face : faces) {
is_iso_face.emplace_back(false);
for (size_t i = 0; i < faces.size(); ++i) {
const auto& face = faces[i];
if (face.supporting_plane > 3) { // plane 0,1,2,3 are tet boundaries
is_iso_face.back() = true;
for (const auto& vid : face.vertices) { is_iso_vert[vid] = true; }
tet_iso_faces.emplace(i);
for (const auto& vId : face.vertices) tet_iso_verts.emplace(vId);
}
}
} else {
for (const auto& face : faces) {
is_iso_face.emplace_back(false);
auto pid = face.supporting_plane;
auto uid = arrangement.unique_plane_indices[pid];
for (const auto& plane_id : arrangement.unique_planes[uid]) {
if (plane_id > 3) { // plane 0,1,2,3 are tet boundaries
is_iso_face.back() = true;
for (const auto& vid : face.vertices) { is_iso_vert[vid] = true; }
break;
for (size_t i = 0; i < faces.size(); ++i) {
const auto& face = faces[i];
const auto pId = face.supporting_plane;
const auto uId = arrangement.unique_plane_indices[pId];
for (const auto& plane_index : arrangement.unique_planes[uId]) {
if (plane_index > 3) { // plane 0,1,2,3 are tet boundaries
tet_iso_faces.emplace(i);
for (const auto& vId : face.vertices) tet_iso_verts.emplace(vId);
}
}
}
}
// map: local vert index --> iso-vert index
iso_vId_of_vert.clear();
iso_vId_of_vert.reserve(vertices.size());
// create iso-vertices
for (uint32_t j = 0; j < vertices.size(); j++) {
iso_vId_of_vert.emplace_back(invalid_index);
if (is_iso_vert[j]) {
uint32_t num_boundary_planes = 0;
uint32_t num_impl_planes = 0;
for (const auto& vertex_plane : vertices[j]) {
if (vertex_plane >= 4) { // plane 0,1,2,3 are tet boundaries
implicit_pIds[num_impl_planes++] =
static_cast<uint16_t>(active_subface_indices[vertex_plane - 4 + start_index]);
} else {
boundary_pIds[num_boundary_planes++] = vertex_plane;
}
for (const auto& vertex_index : tet_iso_verts) {
const auto& vertex = vertices[vertex_index];
auto& local_vert_to_iso_vert = iso_vId_of_vert.emplace_back(invalid_index);
auto& iso_vert = iso_verts.emplace_back();
uint32_t num_boundary_planes{};
uint32_t num_impl_planes{};
for (const auto& vertex_plane : vertex) {
if (vertex_plane >= 4) { // plane 0,1,2,3 are tet boundaries
implicit_pIds[num_impl_planes++] =
static_cast<uint16_t>(active_subface_indices[vertex_plane - 4 + start_index]);
} else {
boundary_pIds[num_boundary_planes++] = vertex_plane;
}
switch (num_boundary_planes) {
case 2: // on tet edge
{
used_pId.fill(false);
used_pId[boundary_pIds[0]] = true;
used_pId[boundary_pIds[1]] = true;
uint32_t num_vIds = 0;
for (uint32_t k = 0; k < 4; k++) {
if (!used_pId[k]) vIds2[num_vIds++] = tetrahedrons[i][k];
}
uint32_t vId1 = vIds2[0];
uint32_t vId2 = vIds2[1];
if (vId1 > vId2) std::swap(vId1, vId2);
key3 = {vId1, vId2, implicit_pIds[0]};
auto iter_inserted = vert_on_tetEdge.try_emplace(key3, static_cast<uint32_t>(iso_verts.size()));
if (iter_inserted.second) {
auto& iso_vert = iso_verts.emplace_back();
iso_vert.header = {i, j, 2};
iso_vert.simplex_vertex_indices = {vId1, vId2};
iso_vert.subface_indices = {implicit_pIds[0]};
}
std::sort(boundary_pIds.begin(), boundary_pIds.begin() + num_boundary_planes);
std::set_difference(tet_boundary_faces.begin(),
tet_boundary_faces.end(),
boundary_pIds.begin(),
boundary_pIds.begin() + num_impl_planes,
not_boundary_pIds.begin());
const auto f1 = vertex_infos(vertex_indices_mapping.at(vId1), implicit_pIds[0]);
const auto f2 = vertex_infos(vertex_indices_mapping.at(vId2), implicit_pIds[0]);
const auto coord = compute_barycentric_coords(f1, f2);
iso_pts.emplace_back(coord[0] * get_vert_pos(vId1, scene_bg_mesh_info)
+ coord[1] * get_vert_pos(vId2, scene_bg_mesh_info));
}
iso_vId_of_vert.back() = iter_inserted.first->second;
break;
}
case 1: // on tet face
{
uint32_t pId = boundary_pIds[0];
uint32_t num_vIds = 0;
for (uint32_t k = 0; k < 4; k++) {
if (k != pId) vIds3[num_vIds++] = tetrahedrons[i][k];
}
std::sort(vIds3.begin(), vIds3.end());
key5 = {vIds3[0], vIds3[1], vIds3[2], implicit_pIds[0], implicit_pIds[1]};
auto iter_inserted = vert_on_tetFace.try_emplace(key5, static_cast<uint32_t>(iso_verts.size()));
if (iter_inserted.second) {
auto& iso_vert = iso_verts.emplace_back();
iso_vert.header = {i, j, 3};
iso_vert.simplex_vertex_indices = {vIds3[0], vIds3[1], vIds3[2]};
iso_vert.subface_indices = {implicit_pIds[0], implicit_pIds[1]};
const auto& tet_vertices = tetrahedrons[tet_index];
iso_vert.header = {tet_index, vertex_index, 4 - num_boundary_planes};
if (num_boundary_planes == 0) { // in tet cell
local_vert_to_iso_vert = static_cast<uint32_t>(iso_verts.size());
iso_vert.simplex_vertex_indices = tet_vertices;
iso_vert.subface_indices = implicit_pIds;
//
const auto &vert_face0 = vertex_infos.row(vertex_indices_mapping.at(vIds3[0])),
vert_face1 = vertex_infos.row(vertex_indices_mapping.at(vIds3[1])),
vert_face2 = vertex_infos.row(vertex_indices_mapping.at(vIds3[2]));
const auto f1 = std::array{vert_face0(implicit_pIds[0]),
vert_face1(implicit_pIds[0]),
vert_face2(implicit_pIds[0])};
const auto f2 = std::array{vert_face0(implicit_pIds[1]),
vert_face1(implicit_pIds[1]),
vert_face2(implicit_pIds[1])};
const auto coord = compute_barycentric_coords(f1, f2);
iso_pts.emplace_back(coord[0] * get_vert_pos(vIds3[0], scene_bg_mesh_info)
+ coord[1] * get_vert_pos(vIds3[1], scene_bg_mesh_info)
+ coord[2] * get_vert_pos(vIds3[2], scene_bg_mesh_info));
}
iso_vId_of_vert.back() = iter_inserted.first->second;
break;
}
case 0: // in tet cell
{
iso_vId_of_vert.back() = static_cast<uint32_t>(iso_verts.size());
auto& iso_vert = iso_verts.emplace_back();
iso_vert.header = {i, j, 4};
iso_vert.simplex_vertex_indices = tetrahedrons[i];
iso_vert.subface_indices = implicit_pIds;
const auto &vert_face0 = vertex_infos.row(vertex_indices_mapping.at(tet_vertices[0])),
vert_face1 = vertex_infos.row(vertex_indices_mapping.at(tet_vertices[1])),
vert_face2 = vertex_infos.row(vertex_indices_mapping.at(tet_vertices[2])),
vert_face3 = vertex_infos.row(vertex_indices_mapping.at(tet_vertices[3]));
const auto f1 = std::array{vert_face0(implicit_pIds[0]),
vert_face1(implicit_pIds[0]),
vert_face2(implicit_pIds[0]),
vert_face3(implicit_pIds[0])};
const auto f2 = std::array{vert_face0(implicit_pIds[1]),
vert_face1(implicit_pIds[1]),
vert_face2(implicit_pIds[1]),
vert_face3(implicit_pIds[1])};
const auto f3 = std::array{vert_face0(implicit_pIds[2]),
vert_face1(implicit_pIds[2]),
vert_face2(implicit_pIds[2]),
vert_face3(implicit_pIds[2])};
const auto coord = compute_barycentric_coords(f1, f2, f3);
iso_pts.emplace_back(coord[0] * get_vert_pos(tet_vertices[0], scene_bg_mesh_info)
+ coord[1] * get_vert_pos(tet_vertices[1], scene_bg_mesh_info)
+ coord[2] * get_vert_pos(tet_vertices[2], scene_bg_mesh_info)
+ coord[3] * get_vert_pos(tet_vertices[3], scene_bg_mesh_info));
} else if (num_boundary_planes == 1) { // on tet face
key5 = {tet_vertices[not_boundary_pIds[0]],
tet_vertices[not_boundary_pIds[1]],
tet_vertices[not_boundary_pIds[2]],
implicit_pIds[0],
implicit_pIds[1]};
auto iter_inserted = vert_on_tetFace.try_emplace(key5, static_cast<uint32_t>(iso_verts.size()));
if (iter_inserted.second) {
iso_vert.simplex_vertex_indices = {key5[0], key5[1], key5[2]};
iso_vert.subface_indices = {implicit_pIds[0], implicit_pIds[1]};
const auto &vert_face0 = vertex_infos.row(vertex_indices_mapping.at(tetrahedrons[i][0])),
vert_face1 = vertex_infos.row(vertex_indices_mapping.at(tetrahedrons[i][1])),
vert_face2 = vertex_infos.row(vertex_indices_mapping.at(tetrahedrons[i][2])),
vert_face3 = vertex_infos.row(vertex_indices_mapping.at(tetrahedrons[i][3]));
const auto f1 = std::array{vert_face0(implicit_pIds[0]),
vert_face1(implicit_pIds[0]),
vert_face2(implicit_pIds[0]),
vert_face3(implicit_pIds[0])};
const auto f2 = std::array{vert_face0(implicit_pIds[1]),
vert_face1(implicit_pIds[1]),
vert_face2(implicit_pIds[1]),
vert_face3(implicit_pIds[1])};
const auto f3 = std::array{vert_face0(implicit_pIds[2]),
vert_face1(implicit_pIds[2]),
vert_face2(implicit_pIds[2]),
vert_face3(implicit_pIds[2])};
const auto coord = compute_barycentric_coords(f1, f2, f3);
iso_pts.emplace_back(coord[0] * get_vert_pos(tetrahedrons[i][0], scene_bg_mesh_info)
+ coord[1] * get_vert_pos(tetrahedrons[i][1], scene_bg_mesh_info)
+ coord[2] * get_vert_pos(tetrahedrons[i][2], scene_bg_mesh_info)
+ coord[3] * get_vert_pos(tetrahedrons[i][3], scene_bg_mesh_info));
//
const auto &vert_face0 = vertex_infos.row(vertex_indices_mapping.at(key5[0])),
vert_face1 = vertex_infos.row(vertex_indices_mapping.at(key5[1])),
vert_face2 = vertex_infos.row(vertex_indices_mapping.at(key5[2]));
const auto f1 =
std::array{vert_face0(implicit_pIds[0]), vert_face1(implicit_pIds[0]), vert_face2(implicit_pIds[0])};
const auto f2 =
std::array{vert_face0(implicit_pIds[1]), vert_face1(implicit_pIds[1]), vert_face2(implicit_pIds[1])};
const auto coord = compute_barycentric_coords(f1, f2);
iso_pts.emplace_back(coord[0] * get_vert_pos(key5[0], scene_bg_mesh_info)
+ coord[1] * get_vert_pos(key5[1], scene_bg_mesh_info)
+ coord[2] * get_vert_pos(key5[2], scene_bg_mesh_info));
}
local_vert_to_iso_vert = iter_inserted.first->second;
} else if (num_boundary_planes == 2) { // on tet edge
key3 = {tet_vertices[not_boundary_pIds[0]], tet_vertices[not_boundary_pIds[1]], implicit_pIds[0]};
auto iter_inserted = vert_on_tetEdge.try_emplace(key3, static_cast<uint32_t>(iso_verts.size()));
if (iter_inserted.second) {
iso_vert.simplex_vertex_indices = {key3[0], key3[1]};
iso_vert.subface_indices = {implicit_pIds[0]};
break;
}
case 3: // on tet vertex
{
used_pId.fill(false);
for (const auto& pId : boundary_pIds) { used_pId[pId] = true; }
uint32_t vId;
for (uint32_t k = 0; k < 4; k++) {
if (!used_pId[k]) {
vId = k;
break;
}
}
auto key = tetrahedrons[i][vId];
auto iter_inserted = vert_on_tetVert.try_emplace(key, static_cast<uint32_t>(iso_verts.size()));
if (iter_inserted.second) {
auto& iso_vert = iso_verts.emplace_back();
iso_vert.header = {i, j, 1};
iso_vert.subface_indices = {static_cast<uint16_t>(tetrahedrons[i][vId])};
const auto f1 = vertex_infos(vertex_indices_mapping.at(key3[0]), implicit_pIds[0]);
const auto f2 = vertex_infos(vertex_indices_mapping.at(key3[1]), implicit_pIds[0]);
const auto coord = compute_barycentric_coords(f1, f2);
iso_pts.emplace_back(coord[0] * get_vert_pos(key3[0], scene_bg_mesh_info)
+ coord[1] * get_vert_pos(key3[1], scene_bg_mesh_info));
}
local_vert_to_iso_vert = iter_inserted.first->second;
} else if (num_boundary_planes == 3) { // on tet vertex
auto key = tet_vertices[not_boundary_pIds[0]];
auto iter_inserted = vert_on_tetVert.try_emplace(key, static_cast<uint32_t>(iso_verts.size()));
if (iter_inserted.second) {
iso_vert.subface_indices = {static_cast<uint16_t>(key)};
iso_pts.emplace_back(get_vert_pos(iso_vert.simplex_vertex_indices[0], scene_bg_mesh_info));
}
iso_vId_of_vert.back() = iter_inserted.first->second;
break;
}
default: break;
iso_pts.emplace_back(get_vert_pos(iso_vert.simplex_vertex_indices[0], scene_bg_mesh_info));
}
local_vert_to_iso_vert = iter_inserted.first->second;
}
}
// create iso-faces
for (uint32_t j = 0; j < faces.size(); j++) {
if (is_iso_face[j]) {
face_verts.clear();
for (unsigned long vId : faces[j].vertices) { face_verts.emplace_back(iso_vId_of_vert[vId]); }
//
// face is on tet boundary if face.negative_cell is NONE
bool face_on_tet_boundary = (faces[j].negative_cell == invalid_index);
//
if (face_on_tet_boundary) {
compute_iso_face_key(face_verts, key3);
auto iter_inserted = face_on_tetFace.try_emplace(key3, static_cast<uint32_t>(iso_faces.size()));
if (iter_inserted.second) {
auto& face = iso_faces.emplace_back();
face.vertex_indices = face_verts;
face.headers.emplace_back(face_header_t{i, j});
face.subface_index = active_subface_indices[faces[j].supporting_plane - 4 + start_index];
} else { // iso_face inserted before
uint32_t iso_face_id = (iter_inserted.first)->second;
iso_faces[iso_face_id].headers.emplace_back(face_header_t{i, j});
}
} else { // face not on tet boundary
auto& face = iso_faces.emplace_back();
face.vertex_indices = face_verts;
face.headers.emplace_back(face_header_t{i, j});
face.subface_index = active_subface_indices[faces[j].supporting_plane - 4 + start_index];
for (const auto& face_index : tet_iso_faces) {
const auto& face = faces[face_index];
face_verts.clear();
for (unsigned long vId : face.vertices) { face_verts.emplace_back(iso_vId_of_vert[vId]); }
const face_header_t face_header{tet_index, face_index};
// face is on tet boundary if face.negative_cell is NONE
uint32_t iso_face_index{static_cast<uint32_t>(iso_faces.size())};
bool is_new_iso_face{};
if (face.negative_cell == invalid_index) {
compute_iso_face_key(face_verts, key3);
auto iter_inserted = face_on_tetFace.try_emplace(key3, static_cast<uint32_t>(iso_faces.size()));
if (iter_inserted.second) {
iso_face_index = (iter_inserted.first)->second;
is_new_iso_face = true;
}
}
if (is_new_iso_face) {
auto& iso_face = iso_faces.emplace_back();
iso_face.vertex_indices = face_verts;
iso_face.subface_index = active_subface_indices[face.supporting_plane - 4 + start_index];
}
iso_faces[iso_face_index].headers.emplace_back(face_header);
}
}
//
}

10
network_process/src/process.cpp
View File

@ -52,6 +52,13 @@ ISNP_API void build_implicit_network_by_blobtree(const s_settings&
tetrahedron_active_subface_start_index,
active_subface_indices,
zero_vertex_to_incident_tet_mapping);
for (auto& [vertex_index, tet_indices] : zero_vertex_to_incident_tet_mapping)
{
std::sort(tet_indices.begin(), tet_indices.end());
auto iter = std::unique(tet_indices.begin(), tet_indices.end());
tet_indices.erase(iter, tet_indices.end());
tet_indices.shrink_to_fit();
}
}
auto tet_active_subface_counts = compute_implicit_arrangements(vertex_infos,
@ -85,7 +92,7 @@ ISNP_API void build_implicit_network_by_blobtree(const s_settings&
stl_vector_mp<uint32_t> shell_to_cell{};
stl_vector_mp<boundary_edge_header_t> chain_edge_headers{};
{
stl_vector_mp<stl_vector_mp<uint32_t>> half_patch_adj_list(2 * patches.size());
stl_vector_mp<stl_vector_mp<uint32_t>> half_patch_adj_list{};
{
stl_vector_mp<chain_header_t> chain_headers{};
{
@ -102,6 +109,7 @@ ISNP_API void build_implicit_network_by_blobtree(const s_settings&
}
// should we take unique part of chain_of_patch?
}
half_patch_adj_list.resize(2 * patches.size());
for (size_t i = 0; i < chain_headers.size(); ++i)
compute_patch_order(tetrahedrons,
chain_headers[i],

Write Preview
Loading…
Cancel
Save