API tools faq
paste
Guest User

mt_implementations

a guest
Mar 5th, 2018
86
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 16.19 KB | None | 0 0
raw download clone embed print report
  1. #include "marching_tets.h"
  2.  
  3. #include <igl/remove_duplicates.h>
  4.  
  5. #include <unordered_map>
  6. #include <vector>
  7. #include <utility>
  8. #include <cstdint>
  9. #include <iostream>
  10. #include <algorithm>
  11.  
  12.  
  13. inline std::int64_t make_edge_key(std::int32_t a, std::int32_t b) {
  14.     std::int64_t ret = 0;
  15.     ret |= a;
  16.     ret |= static_cast<std::int64_t>(b) << 32;
  17.     return ret;
  18. }
  19.  
  20. inline std::int64_t make_edge_key(const std::pair<std::int32_t, std::int32_t>& p) {
  21.     std::int64_t ret = 0;
  22.     ret |= p.first;
  23.     ret |= static_cast<std::int64_t>(p.second) << 32;
  24.     return ret;
  25. }
  26.  
  27. /*
  28.  * Sort based dedup
  29.  */
  30. template <typename DerivedTV,
  31.           typename DerivedTT,
  32.           typename Derivedisovalues,
  33.           typename DerivedoutV,
  34.           typename DerivedoutF>
  35. void igl::marching_tets(
  36.         const Eigen::PlainObjectBase<DerivedTV>& TV,
  37.         const Eigen::PlainObjectBase<DerivedTT>& TT,
  38.         const Eigen::PlainObjectBase<Derivedisovalues>& isovals,
  39.         double isovalue,
  40.         Eigen::PlainObjectBase<DerivedoutV>& outV,
  41.         Eigen::PlainObjectBase<DerivedoutF>& outF) {
  42.     using namespace std;
  43.  
  44.     const int mt_cell_lookup[16][4] = {
  45.         { -1, -1, -1, -1 }, // 0000 0a
  46.         {  0,  2,  1, -1 }, // 0001 1a
  47.         {  0,  3,  4, -1 }, // 0010 1b
  48.         {  2,  1,  3,  4 }, // 0011 2b
  49.         {  5,  3,  1, -1 }, // 0100 1d
  50.         {  0,  2,  5,  3 }, // 0101 2a
  51.         {  0,  1,  5,  4 }, // 0110 2d
  52.         {  2,  5,  4, -1 }, // 0111 3a
  53.         {  4,  5,  2, -1 }, // 1000 1c
  54.         {  0,  4,  5,  1 }, // 1001 2c
  55.         {  0,  3,  5,  2 }, // 1010 2f
  56.         {  1,  3,  5, -1 }, // 1011 3d
  57.         {  4,  3,  1,  2 }, // 1100 2e
  58.         {  0,  4,  3, -1 }, // 1101 3c
  59.         {  0,  1,  2, -1 }, // 1110 3b
  60.         { -1, -1, -1, -1 }, // 1111 0b
  61.     };
  62.  
  63.     const int mt_edge_lookup[6][2] = {
  64.         {0, 1},
  65.         {0, 2},
  66.         {0, 3},
  67.         {1, 2},
  68.         {1, 3},
  69.         {2, 3},
  70.     };
  71.  
  72.     vector<Eigen::RowVector3d> vertices;
  73.     vector<Eigen::RowVector3i> faces;
  74.     vector<pair<pair<int32_t, int32_t>, int32_t>> edge_table;
  75.  
  76.     assert(TT.cols() == 4 && TT.rows() >= 1);
  77.     assert(TV.cols() == 3 && TV.rows() >= 4);
  78.     assert(isovals.cols() == 1);
  79.  
  80.     static_assert(sizeof(std::int64_t) == 2*sizeof(std::int32_t), "need to fit 2 ints into a size_t");
  81.  
  82.     // For each tet
  83.     for (int i = 0; i < TT.rows(); i++) {
  84.         uint8_t key = 0;
  85.         for (int v = 0; v < 4; v++) {
  86.             int vid = TT(i, v);
  87.             uint8_t flag = isovals[vid] > isovalue;
  88.             key |= flag << v;
  89.         }
  90.  
  91.         // Insert the vertices if they don't exist
  92.         int v_ids[4] = {-1, -1, -1, -1}; // This will contain the index in TV of each vertex in the tet
  93.         for (int e = 0; mt_cell_lookup[key][e] != -1 && e < 4; e++) {
  94.             const int v1id = mt_edge_lookup[mt_cell_lookup[key][e]][0];
  95.             const int v2id = mt_edge_lookup[mt_cell_lookup[key][e]][1];
  96.             const int v1i = TT(i, v1id), v2i = TT(i, v2id);
  97.             const Eigen::RowVector3d v1 = TV.row(v1i);
  98.             const Eigen::RowVector3d v2 = TV.row(v2i);
  99.             const double a = fabs(isovals[v1i] - isovalue), b = fabs(isovals[v2i] - isovalue);
  100.             const double w = a / (a+b);
  101.  
  102.             const int vertex_id = vertices.size();
  103.             vertices.push_back((1-w)*v1 + w*v2); // Push back the midpoint
  104.             if (v1i < v2i) {
  105.                 edge_table.push_back(make_pair(make_pair(v1i, v2i), vertex_id));
  106.             } else {
  107.                 edge_table.push_back(make_pair(make_pair(v2i, v1i), vertex_id));
  108.             }
  109.  
  110.             v_ids[e] = vertex_id;
  111.         }
  112.  
  113.         if (v_ids[0] != -1) {
  114.             bool is_quad = mt_cell_lookup[key][3] != -1;
  115.             if (is_quad) {
  116.                 Eigen::RowVector3i f1(v_ids[0], v_ids[1], v_ids[3]);
  117.                 Eigen::RowVector3i f2(v_ids[1], v_ids[2], v_ids[3]);
  118.                 faces.push_back(f1);
  119.                 faces.push_back(f2);
  120.             } else {
  121.                 Eigen::RowVector3i f(v_ids[0], v_ids[1], v_ids[2]);
  122.                 faces.push_back(f);
  123.             }
  124.         }
  125.     }
  126.  
  127.     if(vertices.size() == 0) {
  128.         return;
  129.     }
  130.  
  131.     // Dedup with sort
  132.     std::vector<int> vmap(vertices.size());
  133.  
  134.     outV.resize(vertices.size(), 3);
  135.  
  136.     std::sort(edge_table.begin(), edge_table.end());
  137.     int unique_count = 0;
  138.     for (int i = 0; i < edge_table.size();) { // For each edge
  139.         const pair<int, int> edge = edge_table[i].first;
  140.         outV.row(unique_count) = vertices[edge_table[i].second];
  141.         while (edge == edge_table[i].first) {
  142.             vmap[edge_table[i].second] = unique_count;
  143.             i += 1;
  144.         }
  145.         unique_count += 1;
  146.     }
  147.     outV.conservativeResize(unique_count, 3);
  148.     outF.resize(faces.size(), 3);
  149.     for (int i = 0; i < faces.size(); i++) {
  150.         for (int v = 0; v < 3; v++) {
  151.             outF(i, v) = vmap[faces[i][v]];
  152.         }
  153.     }
  154. }
  155.  
  156. /*
  157.  * Original code
  158.  */
  159. template <typename DerivedTV,
  160.           typename DerivedTT,
  161.           typename Derivedisovalues,
  162.           typename DerivedoutV,
  163.           typename DerivedoutF>
  164. void igl::marching_tets_2(
  165.         const Eigen::PlainObjectBase<DerivedTV>& TV,
  166.         const Eigen::PlainObjectBase<DerivedTT>& TT,
  167.         const Eigen::PlainObjectBase<Derivedisovalues>& isovals,
  168.         double isovalue,
  169.         Eigen::PlainObjectBase<DerivedoutV>& outV,
  170.         Eigen::PlainObjectBase<DerivedoutF>& outF) {
  171.     using namespace std;
  172.  
  173.     const int mt_cell_lookup[16][4] = {
  174.         { -1, -1, -1, -1 }, // 0000 0a
  175.         {  0,  2,  1, -1 }, // 0001 1a
  176.         {  0,  3,  4, -1 }, // 0010 1b
  177.         {  2,  1,  3,  4 }, // 0011 2b
  178.         {  5,  3,  1, -1 }, // 0100 1d
  179.         {  0,  2,  5,  3 }, // 0101 2a
  180.         {  0,  1,  5,  4 }, // 0110 2d
  181.         {  2,  5,  4, -1 }, // 0111 3a
  182.         {  4,  5,  2, -1 }, // 1000 1c
  183.         {  0,  4,  5,  1 }, // 1001 2c
  184.         {  0,  3,  5,  2 }, // 1010 2f
  185.         {  1,  3,  5, -1 }, // 1011 3d
  186.         {  4,  3,  1,  2 }, // 1100 2e
  187.         {  0,  4,  3, -1 }, // 1101 3c
  188.         {  0,  1,  2, -1 }, // 1110 3b
  189.         { -1, -1, -1, -1 }, // 1111 0b
  190.     };
  191.  
  192.     const int mt_edge_lookup[6][2] = {
  193.         {0, 1},
  194.         {0, 2},
  195.         {0, 3},
  196.         {1, 2},
  197.         {1, 3},
  198.         {2, 3},
  199.     };
  200.  
  201.     vector<Eigen::RowVector3d> vertices;
  202.     vector<Eigen::RowVector3i> faces;
  203.     unordered_map<std::int64_t, int> edge_table;
  204.  
  205.     assert(TT.rows() == 4);
  206.     assert(TV.rows() == 3);
  207.     assert(isovals.rows() == 1);
  208.  
  209.     static_assert(sizeof(std::int64_t) == 2*sizeof(std::int32_t), "need to fit 2 ints into a size_t");
  210.  
  211.     // For each tet
  212.     for (int i = 0; i < TT.rows(); i++) {
  213.         uint8_t key = 0;
  214.         for (int v = 0; v < 4; v++) {
  215.             int vid = TT(i, v);
  216.             uint8_t flag = isovals[vid] > isovalue;
  217.             key |= flag << v;
  218.         }
  219.  
  220.         // Insert the vertices if they don't exist
  221.         int v_ids[4] = {-1, -1, -1, -1}; // This will contain the index in TV of each vertex in the tet
  222.         for (int e = 0; mt_cell_lookup[key][e] != -1 && e < 4; e++) {
  223.             std::int64_t edge_table_key = make_edge_key(i, mt_cell_lookup[key][e]);
  224.             auto edge_id_it = edge_table.find(edge_table_key);
  225.  
  226.             if (edge_id_it == edge_table.end()) {
  227.                 const int v1id = mt_edge_lookup[mt_cell_lookup[key][e]][0];
  228.                 const int v2id = mt_edge_lookup[mt_cell_lookup[key][e]][1];
  229.                 const int v1i = TT(i, v1id), v2i = TT(i, v2id);
  230.                 const Eigen::RowVector3d v1 = TV.row(v1i);
  231.                 const Eigen::RowVector3d v2 = TV.row(v2i);
  232.                 const double a = fabs(isovals[v1i] - isovalue), b = fabs(isovals[v2i] - isovalue);
  233.                 const double w = a / (a+b);
  234.                 vertices.push_back((1-w)*v1 + w*v2); // Push back the midpoint
  235.                 edge_table.insert(make_pair(edge_table_key, vertices.size()-1));
  236.                 v_ids[e] = vertices.size()-1;
  237.             } else {
  238.                 v_ids[e] = edge_id_it->second;
  239.             }
  240.         }
  241.         if (v_ids[0] != -1) {
  242.             bool is_quad = mt_cell_lookup[key][3] != -1;
  243.             if (is_quad) {
  244.                 Eigen::RowVector3i f1(v_ids[0], v_ids[1], v_ids[3]);
  245.                 Eigen::RowVector3i f2(v_ids[1], v_ids[2], v_ids[3]);
  246.                 faces.push_back(f1);
  247.                 faces.push_back(f2);
  248.             } else {
  249.                 Eigen::RowVector3i f(v_ids[0], v_ids[1], v_ids[2]);
  250.                 faces.push_back(f);
  251.             }
  252.         }
  253.     }
  254.  
  255.     outV.resize(vertices.size(), 3);
  256.     outF.resize(faces.size(), 3);
  257.     for (int i = 0; i < vertices.size(); i++) {
  258.         outV.row(i) = vertices[i];
  259.     }
  260.     for (int i = 0; i < faces.size(); i++) {
  261.         outF.row(i) = faces[i];
  262.     }
  263. }
  264.  
  265. /*
  266.  * Hash based dedup
  267.  */
  268. template <typename DerivedTV,
  269.           typename DerivedTT,
  270.           typename Derivedisovalues,
  271.           typename DerivedoutV,
  272.           typename DerivedoutF>
  273. void igl::marching_tets_3(
  274.         const Eigen::PlainObjectBase<DerivedTV>& TV,
  275.         const Eigen::PlainObjectBase<DerivedTT>& TT,
  276.         const Eigen::PlainObjectBase<Derivedisovalues>& isovals,
  277.         double isovalue,
  278.         Eigen::PlainObjectBase<DerivedoutV>& outV,
  279.         Eigen::PlainObjectBase<DerivedoutF>& outF) {
  280.     using namespace std;
  281.  
  282.     const int mt_cell_lookup[16][4] = {
  283.         { -1, -1, -1, -1 }, // 0000 0a
  284.         {  0,  2,  1, -1 }, // 0001 1a
  285.         {  0,  3,  4, -1 }, // 0010 1b
  286.         {  2,  1,  3,  4 }, // 0011 2b
  287.         {  5,  3,  1, -1 }, // 0100 1d
  288.         {  0,  2,  5,  3 }, // 0101 2a
  289.         {  0,  1,  5,  4 }, // 0110 2d
  290.         {  2,  5,  4, -1 }, // 0111 3a
  291.         {  4,  5,  2, -1 }, // 1000 1c
  292.         {  0,  4,  5,  1 }, // 1001 2c
  293.         {  0,  3,  5,  2 }, // 1010 2f
  294.         {  1,  3,  5, -1 }, // 1011 3d
  295.         {  4,  3,  1,  2 }, // 1100 2e
  296.         {  0,  4,  3, -1 }, // 1101 3c
  297.         {  0,  1,  2, -1 }, // 1110 3b
  298.         { -1, -1, -1, -1 }, // 1111 0b
  299.     };
  300.  
  301.     const int mt_edge_lookup[6][2] = {
  302.         {0, 1},
  303.         {0, 2},
  304.         {0, 3},
  305.         {1, 2},
  306.         {1, 3},
  307.         {2, 3},
  308.     };
  309.  
  310.     vector<Eigen::RowVector3d> vertices;
  311.     vector<Eigen::RowVector3i> faces;
  312.     vector<pair<int, int>> edge_table;
  313.  
  314.  
  315.     assert(TT.cols() == 4 && TT.rows() >= 1);
  316.     assert(TV.cols() == 3 && TV.rows() >= 4);
  317.     assert(isovals.cols() == 1);
  318.  
  319.     static_assert(sizeof(std::int64_t) == 2*sizeof(std::int32_t), "need to fit 2 ints into a size_t");
  320.  
  321.     // For each tet
  322.     for (int i = 0; i < TT.rows(); i++) {
  323.         uint8_t key = 0;
  324.         for (int v = 0; v < 4; v++) {
  325.             int vid = TT(i, v);
  326.             uint8_t flag = isovals[vid] > isovalue;
  327.             key |= flag << v;
  328.         }
  329.  
  330.         // Insert the vertices if they don't exist
  331.         int v_ids[4] = {-1, -1, -1, -1}; // This will contain the index in TV of each vertex in the tet
  332.         for (int e = 0; mt_cell_lookup[key][e] != -1 && e < 4; e++) {
  333.             const int v1id = mt_edge_lookup[mt_cell_lookup[key][e]][0];
  334.             const int v2id = mt_edge_lookup[mt_cell_lookup[key][e]][1];
  335.             const int v1i = TT(i, v1id), v2i = TT(i, v2id);
  336.             const Eigen::RowVector3d v1 = TV.row(v1i);
  337.             const Eigen::RowVector3d v2 = TV.row(v2i);
  338.             const double a = fabs(isovals[v1i] - isovalue), b = fabs(isovals[v2i] - isovalue);
  339.             const double w = a / (a+b);
  340.  
  341.             const int vertex_id = vertices.size();
  342.             vertices.push_back((1-w)*v1 + w*v2); // Push back the midpoint
  343.             if (v1i < v2i) {
  344.                 edge_table.push_back(make_pair(v1i, v2i));
  345.             } else {
  346.                 edge_table.push_back(make_pair(v2i, v1i));
  347.             }
  348.  
  349.             v_ids[e] = vertex_id;
  350.         }
  351.  
  352.         if (v_ids[0] != -1) {
  353.             bool is_quad = mt_cell_lookup[key][3] != -1;
  354.             if (is_quad) {
  355.                 Eigen::RowVector3i f1(v_ids[0], v_ids[1], v_ids[3]);
  356.                 Eigen::RowVector3i f2(v_ids[1], v_ids[2], v_ids[3]);
  357.                 faces.push_back(f1);
  358.                 faces.push_back(f2);
  359.             } else {
  360.                 Eigen::RowVector3i f(v_ids[0], v_ids[1], v_ids[2]);
  361.                 faces.push_back(f);
  362.             }
  363.         }
  364.     }
  365.  
  366.     // Deduplicate vertices
  367.     int num_unique = 0;
  368.     outV.resize(vertices.size(), 3);
  369.     unordered_map<int64_t, int> emap;
  370.     for (int i = 0; i < faces.size(); i++) {
  371.         for (int v = 0; v < 3; v++) {
  372.             const int vi = faces[i][v];
  373.             const pair<int32_t, int32_t> edge = edge_table[vi];
  374.             const int64_t key = make_edge_key(edge);
  375.             auto it = emap.find(key);
  376.             if (it == emap.end()) { // New unique vertex, insert it
  377.                 outV.row(num_unique) = vertices[vi];
  378.                 outF(i, v) = num_unique;
  379.                 emap.insert(make_pair(key, num_unique));
  380.                 num_unique += 1;
  381.             } else {
  382.                 outF(i, v) = it->second;
  383.             }
  384.         }
  385.     }
  386.     outV.conservativeResize(num_unique, 3);
  387. }
  388.  
  389. #ifdef IGL_STATIC_LIBRARY
  390. template void igl::marching_tets<typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  391.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>,
  392.                                  typename Eigen::Matrix<double, -1, 1, 0, -1, 1>,
  393.                                  typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  394.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>>(
  395.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&,
  396.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&,
  397.         const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&,
  398.         double isovalue,
  399.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
  400.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>&);
  401. template void igl::marching_tets_2<typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  402.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>,
  403.                                  typename Eigen::Matrix<double, -1, 1, 0, -1, 1>,
  404.                                  typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  405.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>>(
  406.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&,
  407.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&,
  408.         const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&,
  409.         double isovalue,
  410.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
  411.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>&);
  412. template void igl::marching_tets_3<typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  413.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>,
  414.                                  typename Eigen::Matrix<double, -1, 1, 0, -1, 1>,
  415.                                  typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  416.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>>(
  417.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&,
  418.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&,
  419.         const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&,
  420.         double isovalue,
  421.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
  422.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>&);
  423. // template void igl::marching_tets<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
  424. #endif // IGL_STATIC_LIBRARY
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!