API tools faq
paste
Guest User

mt_pastebin_fixed

a guest
Mar 5th, 2018
103
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 16.26 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.     outF.resize(faces.size(), 3);
  136.  
  137.     std::sort(edge_table.begin(), edge_table.end());
  138.     int unique_count = 0;
  139.     for (int i = 0; i < edge_table.size();) { // For each edge
  140.         const pair<int, int> edge = edge_table[i].first;
  141.         outV.row(unique_count) = vertices[edge_table[i].second];
  142.         while (edge == edge_table[i].first) {
  143.             vmap[edge_table[i].second] = unique_count;
  144.             i += 1;
  145.         }
  146.         unique_count += 1;
  147.     }
  148.     outV.conservativeResize(unique_count, 3);
  149.     outF.resize(faces.size(), 3);
  150.     for (int i = 0; i < faces.size(); i++) {
  151.         for (int v = 0; v < 3; v++) {
  152.             outF(i, v) = vmap[faces[i][v]];
  153.         }
  154.     }
  155. }
  156.  
  157. /*
  158.  * Original code
  159.  */
  160. template <typename DerivedTV,
  161.           typename DerivedTT,
  162.           typename Derivedisovalues,
  163.           typename DerivedoutV,
  164.           typename DerivedoutF>
  165. void igl::marching_tets_2(
  166.         const Eigen::PlainObjectBase<DerivedTV>& TV,
  167.         const Eigen::PlainObjectBase<DerivedTT>& TT,
  168.         const Eigen::PlainObjectBase<Derivedisovalues>& isovals,
  169.         double isovalue,
  170.         Eigen::PlainObjectBase<DerivedoutV>& outV,
  171.         Eigen::PlainObjectBase<DerivedoutF>& outF) {
  172.     using namespace std;
  173.  
  174.     const int mt_cell_lookup[16][4] = {
  175.         { -1, -1, -1, -1 }, // 0000 0a
  176.         {  0,  2,  1, -1 }, // 0001 1a
  177.         {  0,  3,  4, -1 }, // 0010 1b
  178.         {  2,  1,  3,  4 }, // 0011 2b
  179.         {  5,  3,  1, -1 }, // 0100 1d
  180.         {  0,  2,  5,  3 }, // 0101 2a
  181.         {  0,  1,  5,  4 }, // 0110 2d
  182.         {  2,  5,  4, -1 }, // 0111 3a
  183.         {  4,  5,  2, -1 }, // 1000 1c
  184.         {  0,  4,  5,  1 }, // 1001 2c
  185.         {  0,  3,  5,  2 }, // 1010 2f
  186.         {  1,  3,  5, -1 }, // 1011 3d
  187.         {  4,  3,  1,  2 }, // 1100 2e
  188.         {  0,  4,  3, -1 }, // 1101 3c
  189.         {  0,  1,  2, -1 }, // 1110 3b
  190.         { -1, -1, -1, -1 }, // 1111 0b
  191.     };
  192.  
  193.     const int mt_edge_lookup[6][2] = {
  194.         {0, 1},
  195.         {0, 2},
  196.         {0, 3},
  197.         {1, 2},
  198.         {1, 3},
  199.         {2, 3},
  200.     };
  201.  
  202.     vector<Eigen::RowVector3d> vertices;
  203.     vector<Eigen::RowVector3i> faces;
  204.     unordered_map<std::int64_t, int> edge_table;
  205.  
  206.     assert(TT.rows() == 4);
  207.     assert(TV.rows() == 3);
  208.     assert(isovals.rows() == 1);
  209.  
  210.     static_assert(sizeof(std::int64_t) == 2*sizeof(std::int32_t), "need to fit 2 ints into a size_t");
  211.  
  212.     // For each tet
  213.     for (int i = 0; i < TT.rows(); i++) {
  214.         uint8_t key = 0;
  215.         for (int v = 0; v < 4; v++) {
  216.             int vid = TT(i, v);
  217.             uint8_t flag = isovals[vid] > isovalue;
  218.             key |= flag << v;
  219.         }
  220.  
  221.         // Insert the vertices if they don't exist
  222.         int v_ids[4] = {-1, -1, -1, -1}; // This will contain the index in TV of each vertex in the tet
  223.         for (int e = 0; mt_cell_lookup[key][e] != -1 && e < 4; e++) {
  224.             std::int64_t edge_table_key = make_edge_key(i, mt_cell_lookup[key][e]);
  225.             auto edge_id_it = edge_table.find(edge_table_key);
  226.  
  227.             if (edge_id_it == edge_table.end()) {
  228.                 const int v1id = mt_edge_lookup[mt_cell_lookup[key][e]][0];
  229.                 const int v2id = mt_edge_lookup[mt_cell_lookup[key][e]][1];
  230.                 const int v1i = TT(i, v1id), v2i = TT(i, v2id);
  231.                 const Eigen::RowVector3d v1 = TV.row(v1i);
  232.                 const Eigen::RowVector3d v2 = TV.row(v2i);
  233.                 const double a = fabs(isovals[v1i] - isovalue), b = fabs(isovals[v2i] - isovalue);
  234.                 const double w = a / (a+b);
  235.                 vertices.push_back((1-w)*v1 + w*v2); // Push back the midpoint
  236.                 edge_table.insert(make_pair(edge_table_key, vertices.size()-1));
  237.                 v_ids[e] = vertices.size()-1;
  238.             } else {
  239.                 v_ids[e] = edge_id_it->second;
  240.             }
  241.         }
  242.         if (v_ids[0] != -1) {
  243.             bool is_quad = mt_cell_lookup[key][3] != -1;
  244.             if (is_quad) {
  245.                 Eigen::RowVector3i f1(v_ids[0], v_ids[1], v_ids[3]);
  246.                 Eigen::RowVector3i f2(v_ids[1], v_ids[2], v_ids[3]);
  247.                 faces.push_back(f1);
  248.                 faces.push_back(f2);
  249.             } else {
  250.                 Eigen::RowVector3i f(v_ids[0], v_ids[1], v_ids[2]);
  251.                 faces.push_back(f);
  252.             }
  253.         }
  254.     }
  255.  
  256.     outV.resize(vertices.size(), 3);
  257.     outF.resize(faces.size(), 3);
  258.     for (int i = 0; i < vertices.size(); i++) {
  259.         outV.row(i) = vertices[i];
  260.     }
  261.     for (int i = 0; i < faces.size(); i++) {
  262.         outF.row(i) = faces[i];
  263.     }
  264. }
  265.  
  266. /*
  267.  * Hash based dedup
  268.  */
  269. template <typename DerivedTV,
  270.           typename DerivedTT,
  271.           typename Derivedisovalues,
  272.           typename DerivedoutV,
  273.           typename DerivedoutF>
  274. void igl::marching_tets_3(
  275.         const Eigen::PlainObjectBase<DerivedTV>& TV,
  276.         const Eigen::PlainObjectBase<DerivedTT>& TT,
  277.         const Eigen::PlainObjectBase<Derivedisovalues>& isovals,
  278.         double isovalue,
  279.         Eigen::PlainObjectBase<DerivedoutV>& outV,
  280.         Eigen::PlainObjectBase<DerivedoutF>& outF) {
  281.     using namespace std;
  282.  
  283.     const int mt_cell_lookup[16][4] = {
  284.         { -1, -1, -1, -1 }, // 0000 0a
  285.         {  0,  2,  1, -1 }, // 0001 1a
  286.         {  0,  3,  4, -1 }, // 0010 1b
  287.         {  2,  1,  3,  4 }, // 0011 2b
  288.         {  5,  3,  1, -1 }, // 0100 1d
  289.         {  0,  2,  5,  3 }, // 0101 2a
  290.         {  0,  1,  5,  4 }, // 0110 2d
  291.         {  2,  5,  4, -1 }, // 0111 3a
  292.         {  4,  5,  2, -1 }, // 1000 1c
  293.         {  0,  4,  5,  1 }, // 1001 2c
  294.         {  0,  3,  5,  2 }, // 1010 2f
  295.         {  1,  3,  5, -1 }, // 1011 3d
  296.         {  4,  3,  1,  2 }, // 1100 2e
  297.         {  0,  4,  3, -1 }, // 1101 3c
  298.         {  0,  1,  2, -1 }, // 1110 3b
  299.         { -1, -1, -1, -1 }, // 1111 0b
  300.     };
  301.  
  302.     const int mt_edge_lookup[6][2] = {
  303.         {0, 1},
  304.         {0, 2},
  305.         {0, 3},
  306.         {1, 2},
  307.         {1, 3},
  308.         {2, 3},
  309.     };
  310.  
  311.     vector<Eigen::RowVector3d> vertices;
  312.     vector<Eigen::RowVector3i> faces;
  313.     vector<pair<int, int>> edge_table;
  314.  
  315.  
  316.     assert(TT.cols() == 4 && TT.rows() >= 1);
  317.     assert(TV.cols() == 3 && TV.rows() >= 4);
  318.     assert(isovals.cols() == 1);
  319.  
  320.     static_assert(sizeof(std::int64_t) == 2*sizeof(std::int32_t), "need to fit 2 ints into a size_t");
  321.  
  322.     // For each tet
  323.     for (int i = 0; i < TT.rows(); i++) {
  324.         uint8_t key = 0;
  325.         for (int v = 0; v < 4; v++) {
  326.             int vid = TT(i, v);
  327.             uint8_t flag = isovals[vid] > isovalue;
  328.             key |= flag << v;
  329.         }
  330.  
  331.         // Insert the vertices if they don't exist
  332.         int v_ids[4] = {-1, -1, -1, -1}; // This will contain the index in TV of each vertex in the tet
  333.         for (int e = 0; mt_cell_lookup[key][e] != -1 && e < 4; e++) {
  334.             const int v1id = mt_edge_lookup[mt_cell_lookup[key][e]][0];
  335.             const int v2id = mt_edge_lookup[mt_cell_lookup[key][e]][1];
  336.             const int v1i = TT(i, v1id), v2i = TT(i, v2id);
  337.             const Eigen::RowVector3d v1 = TV.row(v1i);
  338.             const Eigen::RowVector3d v2 = TV.row(v2i);
  339.             const double a = fabs(isovals[v1i] - isovalue), b = fabs(isovals[v2i] - isovalue);
  340.             const double w = a / (a+b);
  341.  
  342.             const int vertex_id = vertices.size();
  343.             vertices.push_back((1-w)*v1 + w*v2); // Push back the midpoint
  344.             if (v1i < v2i) {
  345.                 edge_table.push_back(make_pair(v1i, v2i));
  346.             } else {
  347.                 edge_table.push_back(make_pair(v2i, v1i));
  348.             }
  349.  
  350.             v_ids[e] = vertex_id;
  351.         }
  352.  
  353.         if (v_ids[0] != -1) {
  354.             bool is_quad = mt_cell_lookup[key][3] != -1;
  355.             if (is_quad) {
  356.                 Eigen::RowVector3i f1(v_ids[0], v_ids[1], v_ids[3]);
  357.                 Eigen::RowVector3i f2(v_ids[1], v_ids[2], v_ids[3]);
  358.                 faces.push_back(f1);
  359.                 faces.push_back(f2);
  360.             } else {
  361.                 Eigen::RowVector3i f(v_ids[0], v_ids[1], v_ids[2]);
  362.                 faces.push_back(f);
  363.             }
  364.         }
  365.     }
  366.  
  367.     // Deduplicate vertices
  368.     int num_unique = 0;
  369.     outV.resize(vertices.size(), 3);
  370.     outF.resize(faces.size(), 3);
  371.     unordered_map<int64_t, int> emap;
  372.     for (int i = 0; i < faces.size(); i++) {
  373.         for (int v = 0; v < 3; v++) {
  374.             const int vi = faces[i][v];
  375.             const pair<int32_t, int32_t> edge = edge_table[vi];
  376.             const int64_t key = make_edge_key(edge);
  377.             auto it = emap.find(key);
  378.             if (it == emap.end()) { // New unique vertex, insert it
  379.                 outV.row(num_unique) = vertices[vi];
  380.                 outF(i, v) = num_unique;
  381.                 emap.insert(make_pair(key, num_unique));
  382.                 num_unique += 1;
  383.             } else {
  384.                 outF(i, v) = it->second;
  385.             }
  386.         }
  387.     }
  388.     outV.conservativeResize(num_unique, 3);
  389. }
  390.  
  391. #ifdef IGL_STATIC_LIBRARY
  392. template void igl::marching_tets<typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  393.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>,
  394.                                  typename Eigen::Matrix<double, -1, 1, 0, -1, 1>,
  395.                                  typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  396.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>>(
  397.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&,
  398.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&,
  399.         const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&,
  400.         double isovalue,
  401.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
  402.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>&);
  403. template void igl::marching_tets_2<typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  404.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>,
  405.                                  typename Eigen::Matrix<double, -1, 1, 0, -1, 1>,
  406.                                  typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  407.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>>(
  408.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&,
  409.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&,
  410.         const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&,
  411.         double isovalue,
  412.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
  413.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>&);
  414. template void igl::marching_tets_3<typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  415.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>,
  416.                                  typename Eigen::Matrix<double, -1, 1, 0, -1, 1>,
  417.                                  typename Eigen::Matrix<double, -1, -1, 0, -1, -1>,
  418.                                  typename Eigen::Matrix<int, -1, -1, 0, -1, -1>>(
  419.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&,
  420.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&,
  421.         const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&,
  422.         double isovalue,
  423.         Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
  424.         Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>&);
  425. // 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> >&);
  426. #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!