| 55 | } |
| 56 | |
| 57 | VectorF VoxelUtils::is_delaunay( |
| 58 | const MatrixFr& vertices, |
| 59 | const MatrixIr& tets) { |
| 60 | const size_t num_tets = tets.rows(); |
| 61 | VectorF result(num_tets); |
| 62 | result.setConstant(1); |
| 63 | |
| 64 | MultipletMap<Triplet, size_t> adj_list; |
| 65 | for (size_t i=0; i<num_tets; i++) { |
| 66 | adj_list.insert({tets(i,0), tets(i,1), tets(i,2)}, i); |
| 67 | adj_list.insert({tets(i,1), tets(i,2), tets(i,3)}, i); |
| 68 | adj_list.insert({tets(i,2), tets(i,3), tets(i,0)}, i); |
| 69 | adj_list.insert({tets(i,3), tets(i,0), tets(i,1)}, i); |
| 70 | } |
| 71 | |
| 72 | constexpr int INVALID = std::numeric_limits<int>::max(); |
| 73 | auto get_opposite_vertex = [&tets](size_t index, const Triplet& f) { |
| 74 | for (size_t i=0; i<4; i++) { |
| 75 | const int val = tets(index, i); |
| 76 | if (val != f.get_data()[0] && |
| 77 | val != f.get_data()[1] && |
| 78 | val != f.get_data()[2]) { |
| 79 | return val; |
| 80 | } |
| 81 | } |
| 82 | // Tet is topologically degenerate. |
| 83 | return INVALID; |
| 84 | }; |
| 85 | |
| 86 | for (size_t i=0; i<num_tets; i++) { |
| 87 | const int i0 = tets(i,0); |
| 88 | const int i1 = tets(i,1); |
| 89 | const int i2 = tets(i,2); |
| 90 | const int i3 = tets(i,3); |
| 91 | const Vector3F v0 = vertices.row(i0); |
| 92 | const Vector3F v1 = vertices.row(i1); |
| 93 | const Vector3F v2 = vertices.row(i2); |
| 94 | const Vector3F v3 = vertices.row(i3); |
| 95 | |
| 96 | std::vector<Triplet> faces = { |
| 97 | {i0, i1, i2}, |
| 98 | {i1, i2, i3}, |
| 99 | {i2, i3, i0}, |
| 100 | {i3, i0, i1} |
| 101 | }; |
| 102 | |
| 103 | for (const auto& f : faces) { |
| 104 | if (result[i] < 0) break; |
| 105 | const auto& adj_tets = adj_list[f]; |
| 106 | for (auto j : adj_tets) { |
| 107 | if (i==j) continue; |
| 108 | if (result[i] < 0) break; |
| 109 | const int oppo = get_opposite_vertex(j, f); |
| 110 | const Vector3F vo = vertices.row(oppo); |
| 111 | // Note that the orientation of the sphere/tet is different from |
| 112 | // the orientation defined by the MSH format. Swapping v0 and |
| 113 | // v1 to ensure consistency. |
| 114 | auto r = insphere( |
no test coverage detected