diff --git a/network_process/interface/fwd_types.hpp b/network_process/interface/fwd_types.hpp index 9e52f7e..11e978c 100644 --- a/network_process/interface/fwd_types.hpp +++ b/network_process/interface/fwd_types.hpp @@ -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(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(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(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); diff --git a/network_process/src/connect_by_topo/pair_faces.cpp b/network_process/src/connect_by_topo/pair_faces.cpp index 2b4b244..c0305aa 100644 --- a/network_process/src/connect_by_topo/pair_faces.cpp +++ b/network_process/src/connect_by_topo/pair_faces.cpp @@ -1,5 +1,12 @@ +#include +#include #include +/// 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> &tetrahedrons, const chain_header_t &character_edge, const stl_vector_mp &iso_verts, @@ -28,42 +35,53 @@ void compute_patch_order(const stl_vector_mp> 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 incident_tets1(tet_ids.begin(), tet_ids.end()); - stl_vector_mp 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 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> } 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 tet1_vIds(tetrahedrons[tet_id1].begin(), tetrahedrons[tet_id1].end()); - stl_vector_mp 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 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> } // =============================================================================================== +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> &half_patch_adj_list) { // find tet faces that are incident to the character_edge - stl_vector_mp 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 iso_face_Id_of_face(tet_cut_result.faces.size(), invalid_index); + stl_vector_mp 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 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 &containing_simplex, @@ -206,38 +222,30 @@ void pair_patches_in_tets(const stl_vector_mp &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 iso_face_Id_of_face{}; + flat_hash_map_mp 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 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 &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 opposite_face{}; - // map: (tet_Id, tet_vert_Id) -> boundary vert Id - flat_hash_map_mp boundary_vert_id{}; - uint32_t num_boundary_vert = 0; - // hash table for vertices on the boundary of tetrahedron - flat_hash_map_mp vert_on_tetVert{}; - flat_hash_map_mp, uint32_t> vert_on_tetEdge{}; - flat_hash_map_mp, uint32_t> vert_on_tetFace{}; + // keys of vertices on the boundary of tetrahedron + stl_vector_mp tet_vert_keys{}; + stl_vector_mp> tet_edge_keys{}; + stl_vector_mp> 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, face_header_t> face_on_tetFace{}; // auxiliary data - stl_vector_mp is_boundary_vert{}; - stl_vector_mp is_boundary_face{}; + flat_hash_set_mp boundary_vertices{}; + flat_hash_set_mp boundary_faces{}; + // map: local vert index --> boundary vert index stl_vector_mp boundary_vId_of_vert{}; stl_vector_mp face_verts{}; pod_key_t<3> key3; pod_key_t<5> key5; - std::array used_pId; - std::array vIds2; - std::array vIds3; std::array implicit_pIds; std::array boundary_pIds; - for (const auto i : containing_tetIds) { + std::array 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::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; } } \ No newline at end of file diff --git a/network_process/src/connect_by_topo/patch_connectivity.cpp b/network_process/src/connect_by_topo/patch_connectivity.cpp index 35a59fc..30357d0 100644 --- a/network_process/src/connect_by_topo/patch_connectivity.cpp +++ b/network_process/src/connect_by_topo/patch_connectivity.cpp @@ -13,7 +13,7 @@ void compute_patch_edges(const stl_vector_mp& patch_face edges_of_face.reserve(patch_faces.size()); uint32_t num_iso_edge{}; // map: (v1, v2) -> iso-edge index - flat_hash_map_mp, uint32_t> edge_id{}; + // flat_hash_map_mp, 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(face.vertex_indices.size()); diff --git a/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp b/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp index 25247b4..9079aeb 100644 --- a/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp +++ b/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp @@ -21,78 +21,139 @@ void build_tetrahedron_and_adjacency(const scene_bg_mesh_info_t& flat_hash_map_mp>& reverse_vertex_adjacency, btree_map_mp>& vertex_to_tet_mapping) { - auto max_pos = scene_bg_mesh_info.vert_resolution.template cast() - Eigen::Vector::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::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); diff --git a/network_process/src/prim_gen/extract_mesh.cpp b/network_process/src/prim_gen/extract_mesh.cpp index 39bd83e..7048f15 100644 --- a/network_process/src/prim_gen/extract_mesh.cpp +++ b/network_process/src/prim_gen/extract_mesh.cpp @@ -1,6 +1,8 @@ #include #include +const std::array tet_boundary_faces = {0u, 1u, 2u, 3u}; + void extract_iso_mesh(const std::array& tet_active_subface_counts, const stl_vector_mp& tetrahedron_arrangements, const scene_bg_mesh_info_t& scene_bg_mesh_info, @@ -32,225 +34,179 @@ void extract_iso_mesh(const std::array& tet_active_s flat_hash_map_mp, uint32_t> face_on_tetFace{}; // - stl_vector_mp is_iso_vert{}; - is_iso_vert.reserve(8); - stl_vector_mp is_iso_face{}; - is_iso_face.reserve(9); - stl_vector_mp iso_vId_of_vert{}; - iso_vId_of_vert.reserve(8); - stl_vector_mp face_verts{}; + flat_hash_set_mp tet_iso_verts{}; + flat_hash_set_mp tet_iso_faces{}; + // map: local vert index --> iso-vert index + stl_vector_mp iso_vId_of_vert{}; + stl_vector_mp face_verts{}; face_verts.reserve(4); pod_key_t<3> key3; pod_key_t<5> key5; - std::array used_pId; - std::array vIds2; - std::array vIds3; std::array implicit_pIds; std::array boundary_pIds; + std::array 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(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(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(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(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(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(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(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(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(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(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(iso_verts.size())); + if (iter_inserted.second) { + iso_vert.subface_indices = {static_cast(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(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(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(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); } } - // } \ No newline at end of file diff --git a/network_process/src/process.cpp b/network_process/src/process.cpp index c7faef5..41c5c51 100644 --- a/network_process/src/process.cpp +++ b/network_process/src/process.cpp @@ -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 shell_to_cell{}; stl_vector_mp chain_edge_headers{}; { - stl_vector_mp> half_patch_adj_list(2 * patches.size()); + stl_vector_mp> half_patch_adj_list{}; { stl_vector_mp 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],