Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- template <typename DerivedTV,
- typename DerivedTT,
- typename Derivedisovalues,
- typename DerivedoutV,
- typename DerivedoutF>
- void igl::marching_tets_4(
- const Eigen::PlainObjectBase<DerivedTV>& TV,
- const Eigen::PlainObjectBase<DerivedTT>& TT,
- const Eigen::PlainObjectBase<Derivedisovalues>& isovals,
- double isovalue,
- Eigen::PlainObjectBase<DerivedoutV>& outV,
- Eigen::PlainObjectBase<DerivedoutF>& outF) {
- const int mt_cell_lookup[16][4] = {
- { -1, -1, -1, -1 }, // 0000 0a
- { 0, 2, 1, -1 }, // 0001 1a
- { 0, 3, 4, -1 }, // 0010 1b
- { 2, 1, 3, 4 }, // 0011 2b
- { 5, 3, 1, -1 }, // 0100 1d
- { 0, 2, 5, 3 }, // 0101 2a
- { 0, 1, 5, 4 }, // 0110 2d
- { 2, 5, 4, -1 }, // 0111 3a
- { 4, 5, 2, -1 }, // 1000 1c
- { 0, 4, 5, 1 }, // 1001 2c
- { 0, 3, 5, 2 }, // 1010 2f
- { 1, 3, 5, -1 }, // 1011 3d
- { 4, 3, 1, 2 }, // 1100 2e
- { 0, 4, 3, -1 }, // 1101 3c
- { 0, 1, 2, -1 }, // 1110 3b
- { -1, -1, -1, -1 }, // 1111 0b
- };
- const int mt_edge_lookup[6][2] = {
- {0, 1},
- {0, 2},
- {0, 3},
- {1, 2},
- {1, 3},
- {2, 3},
- };
- typedef std::pair<int, int> Edge;
- std::vector<Eigen::RowVector3i> faces;
- std::vector<std::pair<Edge, int>> all_edges;
- auto make_edge = [](int v1, int v2) {
- return std::make_pair(std::min(v1, v2), std::max(v1, v2));
- };
- assert(TT.rows() == 4);
- assert(TV.rows() == 3);
- assert(isovals.rows() == 1);
- static_assert(sizeof(std::int64_t) == 2*sizeof(std::int32_t), "need to fit 2 ints into a size_t");
- // For each tet
- for (int i = 0; i < TT.rows(); i++) {
- uint8_t key = 0;
- for (int v = 0; v < 4; v++) {
- int vid = TT(i, v);
- uint8_t flag = isovals[vid] > isovalue;
- key |= flag << v;
- }
- // Insert the vertices if they don't exist
- int v_ids[4] = {-1, -1, -1, -1}; // This will contain the index in TV of each vertex in the tet
- for (int e = 0; e < 4 && mt_cell_lookup[key][e] != -1; ++e) {
- int v1id = mt_edge_lookup[mt_cell_lookup[key][e]][0];
- int v2id = mt_edge_lookup[mt_cell_lookup[key][e]][1];
- int v1 = TT(i, v1id);
- int v2 = TT(i, v2id);
- v_ids[e] = all_edges.size();
- all_edges.emplace_back(make_edge(v1, v2), v_ids[e]);
- }
- if (v_ids[0] != -1) {
- bool is_quad = mt_cell_lookup[key][3] != -1;
- if (is_quad) {
- Eigen::RowVector3i f1(v_ids[0], v_ids[1], v_ids[3]);
- Eigen::RowVector3i f2(v_ids[1], v_ids[2], v_ids[3]);
- faces.push_back(f1);
- faces.push_back(f2);
- } else {
- Eigen::RowVector3i f(v_ids[0], v_ids[1], v_ids[2]);
- faces.push_back(f);
- }
- }
- }
- // Dedup vertices
- Eigen::VectorXi remap(all_edges.size());
- std::sort(all_edges.begin(), all_edges.end());
- Edge prev(-1, -1);
- int num_vertices = -1;
- for (auto &kv : all_edges) {
- if (prev != kv.first) { ++num_vertices; }
- remap(kv.second) = num_vertices;
- prev = kv.first;
- }
- ++num_vertices;
- // Assign final mesh
- outV.resize(num_vertices, 3);
- outF.resize(faces.size(), 3);
- for (auto &kv : all_edges) {
- int v1 = kv.first.first;
- int v2 = kv.first.second;
- Eigen::RowVector3d p1 = TV.row(v1);
- Eigen::RowVector3d p2 = TV.row(v2);
- auto s0 = fabs(isovals[v1] - isovalue);
- auto s1 = fabs(isovals[v2] - isovalue);
- auto t = s0 / (s0+s1);
- outV.row(remap(kv.second)) = ((1.0f-t)*p1 + t*p2);
- }
- for (int i = 0; i < faces.size(); i++) {
- outF.row(i) = faces[i].unaryExpr([&](int v) { return remap(v); });
- }
- }
Advertisement
Add Comment
Please, Sign In to add comment
