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_3(
- 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) {
- using namespace std;
- const auto make_edge_key = [](const std::pair<std::int32_t, std::int32_t>& p) -> std::int64_t {
- static_assert(sizeof(std::int64_t) == 2*sizeof(std::int32_t), "need to fit 2 ints into a size_t");
- std::int64_t ret = 0;
- ret |= p.first;
- ret |= static_cast<std::int64_t>(p.second) << 32;
- return ret;
- };
- 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},
- };
- // vector<Eigen::RowVector3d> vertices;
- vector<Eigen::RowVector3i> faces;
- vector<pair<int, int>> edge_table;
- assert(TT.cols() == 4 && TT.rows() >= 1);
- assert(TV.cols() == 3 && TV.rows() >= 4);
- assert(isovals.cols() == 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++) {
- const int v1id = mt_edge_lookup[mt_cell_lookup[key][e]][0];
- const int v2id = mt_edge_lookup[mt_cell_lookup[key][e]][1];
- const int v1i = TT(i, v1id), v2i = TT(i, v2id);
- const int vertex_id = edge_table.size();
- edge_table.push_back(make_pair(min(v1i, v2i), max(v1i, v2i)));
- v_ids[e] = vertex_id;
- }
- 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);
- }
- }
- }
- // Deduplicate vertices
- int num_unique = 0;
- outV.resize(edge_table.size(), 3);
- outF.resize(faces.size(), 3);
- unordered_map<int64_t, int> emap;
- emap.max_load_factor(0.5);
- emap.reserve(edge_table.size());
- for (int i = 0; i < faces.size(); i++) {
- for (int v = 0; v < 3; v++) {
- const int vi = faces[i][v];
- const pair<int32_t, int32_t> edge = edge_table[vi];
- const int64_t key = make_edge_key(edge);
- auto it = emap.find(key);
- if (it == emap.end()) { // New unique vertex, insert it
- const Eigen::RowVector3d v1 = TV.row(edge.first);
- const Eigen::RowVector3d v2 = TV.row(edge.second);
- const double a = fabs(isovals[edge.first] - isovalue), b = fabs(isovals[edge.second] - isovalue);
- const double w = a / (a+b);
- outV.row(num_unique) = (1-w)*v1 + w*v2;
- outF(i, v) = num_unique;
- emap.insert(make_pair(key, num_unique));
- num_unique += 1;
- } else {
- outF(i, v) = it->second;
- }
- }
- }
- outV.conservativeResize(num_unique, 3);
- }
Advertisement
Add Comment
Please, Sign In to add comment
