Marching Tetrahedra (#716)

* Marching Tetrahedra

* dropping cxx11
This commit is contained in:
Francis Williams
2018-10-14 19:14:45 -04:00
committed by Alec Jacobson
parent b4eaa0dac7
commit d152b6072e
17 changed files with 462 additions and 467 deletions
+194
View File
@@ -0,0 +1,194 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2018 Francis Williams <francis@fwilliams.info>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include "marching_tets.h"
#include <unordered_map>
#include <vector>
#include <utility>
#include <cstdint>
#include <iostream>
template <typename DerivedTV,
typename DerivedTT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ,
typename BCType>
void igl::marching_tets(
const Eigen::PlainObjectBase<DerivedTV>& TV,
const Eigen::PlainObjectBase<DerivedTT>& TT,
const Eigen::PlainObjectBase<DerivedS>& isovals,
double isovalue,
Eigen::PlainObjectBase<DerivedSV>& outV,
Eigen::PlainObjectBase<DerivedSF>& outF,
Eigen::PlainObjectBase<DerivedJ>& J,
Eigen::SparseMatrix<BCType>& BC)
{
using namespace std;
// We're hashing edges to deduplicate using 64 bit ints. The upper and lower
// 32 bits of a key are the indices of vertices in the mesh. The implication is
// that you can only have 2^32 vertices which I have deemed sufficient for
// anything reasonable.
const auto make_edge_key = [](const pair<int32_t, int32_t>& p) -> int64_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 },
{ 0, 2, 1, -1 },
{ 0, 3, 4, -1 },
{ 2, 1, 3, 4 },
{ 5, 3, 1, -1 },
{ 0, 2, 5, 3 },
{ 0, 1, 5, 4 },
{ 2, 5, 4, -1 },
{ 4, 5, 2, -1 },
{ 0, 4, 5, 1 },
{ 0, 3, 5, 2 },
{ 1, 3, 5, -1 },
{ 4, 3, 1, 2 },
{ 0, 4, 3, -1 },
{ 0, 1, 2, -1 },
{ -1, -1, -1, -1 },
};
const int mt_edge_lookup[6][2] =
{
{0, 1},
{0, 2},
{0, 3},
{1, 2},
{1, 3},
{2, 3},
};
// Store the faces and the tet they are in
vector<pair<Eigen::RowVector3i, int>> faces;
// Store the edges in the tet mesh which we add vertices on
// so we can deduplicate
vector<pair<int, int>> edge_table;
assert(TT.cols() == 4 && TT.rows() >= 1);
assert(TV.cols() == 3 && TV.rows() >= 4);
assert(isovals.cols() == 1);
// For each tet
for (int i = 0; i < TT.rows(); i++)
{
uint8_t key = 0;
for (int v = 0; v < 4; v++)
{
const int vid = TT(i, v);
const uint8_t flag = isovals[vid] > isovalue;
key |= flag << v;
}
// This will contain the index in TV of each vertex in the tet
int v_ids[4] = {-1, -1, -1, -1};
// Insert any vertices if the tet intersects the level surface
for (int e = 0; e < 4 && mt_cell_lookup[key][e] != -1; e++)
{
const int tv1_idx = TT(i, mt_edge_lookup[mt_cell_lookup[key][e]][0]);
const int tv2_idx = TT(i, mt_edge_lookup[mt_cell_lookup[key][e]][1]);
const int vertex_id = edge_table.size();
edge_table.push_back(make_pair(min(tv1_idx, tv2_idx), max(tv1_idx, tv2_idx)));
v_ids[e] = vertex_id;
}
// Insert the corresponding faces
if (v_ids[0] != -1)
{
bool is_quad = mt_cell_lookup[key][3] != -1;
if (is_quad)
{
const Eigen::RowVector3i f1(v_ids[0], v_ids[1], v_ids[3]);
const Eigen::RowVector3i f2(v_ids[1], v_ids[2], v_ids[3]);
faces.push_back(make_pair(f1, i));
faces.push_back(make_pair(f2, i));
}
else
{
const Eigen::RowVector3i f(v_ids[0], v_ids[1], v_ids[2]);
faces.push_back(make_pair(f, i));
}
}
}
int num_unique = 0;
outV.resize(edge_table.size(), 3);
outF.resize(faces.size(), 3);
J.resize(faces.size());
// Sparse matrix triplets for BC
vector<Eigen::Triplet<BCType>> bc_triplets;
bc_triplets.reserve(edge_table.size());
// Deduplicate vertices
unordered_map<int64_t, int> emap;
emap.max_load_factor(0.5);
emap.reserve(edge_table.size());
for (int f = 0; f < faces.size(); f++)
{
for (int v = 0; v < 3; v++)
{
const int vi = faces[f].first[v];
const int ti = faces[f].second;
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
{
// Typedef to make sure we handle floats properly
typedef Eigen::Matrix<typename DerivedTV::Scalar, 1, 3, Eigen::RowMajor, 1, 3> RowVector;
const RowVector v1 = TV.row(edge.first);
const RowVector v2 = TV.row(edge.second);
const double a = fabs(isovals[edge.first] - isovalue);
const double b = fabs(isovals[edge.second] - isovalue);
const double w = a / (a+b);
// Create a casted copy in case BCType is a float and we need to downcast
const BCType bc_w = static_cast<BCType>(w);
bc_triplets.push_back(Eigen::Triplet<BCType>(num_unique, edge.first, 1-bc_w));
bc_triplets.push_back(Eigen::Triplet<BCType>(num_unique, edge.second, bc_w));
// Create a casted copy in case DerivedTV::Scalar is a float and we need to downcast
const typename DerivedTV::Scalar v_w = static_cast<typename DerivedTV::Scalar>(w);
outV.row(num_unique) = (1-v_w)*v1 + v_w*v2;
outF(f, v) = num_unique;
J[f] = ti;
emap.emplace(key, num_unique);
num_unique += 1;
} else {
outF(f, v) = it->second;
}
}
}
outV.conservativeResize(num_unique, 3);
J.conservativeResize(num_unique, 1);
BC.resize(num_unique, TV.rows());
BC.setFromTriplets(bc_triplets.begin(), bc_triplets.end());
}
#ifdef IGL_STATIC_LIBRARY
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::Matrix<int, -1, 1, 0, -1, 1>, double>(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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::SparseMatrix<double, 0, int>&);
#endif // IGL_STATIC_LIBRARY
+196
View File
@@ -0,0 +1,196 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2018 Francis Williams <francis@fwilliams.info>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#ifndef IGL_MARCHING_TETS_H
#define IGL_MARCHING_TETS_H
#include "igl_inline.h"
#include <Eigen/Core>
#include <Eigen/Sparse>
namespace igl {
// marching_tets( TV, TT, S, isovalue, SV, SF, J, BC)
//
// performs the marching tetrahedra algorithm on a tet mesh defined by TV and
// TT with scalar values defined at each vertex in TV. The output is a
// triangle mesh approximating the isosurface coresponding to the value
// isovalue.
//
// Input:
// TV #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
// TT #tets x 4 array -- The indexes of each tet in the tetrahedral mesh
// S #tet_vertices x 1 array -- The values defined on each tet vertex
// isovalue scalar -- The isovalue of the level set we want to compute
//
// Output:
// SV #SV x 3 array -- The vertices of the output level surface mesh
// SF #SF x 3 array -- The face indexes of the output level surface mesh
// J #SF list of indices into TT revealing which tet each face comes from
// BC #SV x #TV list of barycentric coordinates so that SV = BC*TV
template <typename DerivedTV,
typename DerivedTT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ,
typename BCType>
IGL_INLINE void marching_tets(
const Eigen::PlainObjectBase<DerivedTV>& TV,
const Eigen::PlainObjectBase<DerivedTT>& TT,
const Eigen::PlainObjectBase<DerivedS>& S,
double isovalue,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J,
Eigen::SparseMatrix<BCType>& BC);
// marching_tets( TV, TT, S, SV, SF, J, BC)
//
// Performs the marching tetrahedra algorithm on a tet mesh defined by TV and
// TT with scalar values defined at each vertex in TV. The output is a
// triangle mesh approximating the isosurface coresponding to an isovalue of 0.
//
// Input:
// TV #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
// TT #tets x 4 array -- The indexes of each tet in the tetrahedral mesh
// S #tet_vertices x 1 array -- The values defined on each tet vertex
// isovalue scalar -- The isovalue of the level set we want to compute
//
// Output:
// SV #SV x 3 array -- The vertices of the output level surface mesh
// SF #SF x 3 array -- The face indexes of the output level surface mesh
// J #SF list of indices into TT revealing which tet each face comes from
// BC #SV x #TV list of barycentric coordinates so that SV = BC*TV
template <typename DerivedTV,
typename DerivedTT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ,
typename BCType>
IGL_INLINE void marching_tets(
const Eigen::PlainObjectBase<DerivedTV>& TV,
const Eigen::PlainObjectBase<DerivedTT>& TT,
const Eigen::PlainObjectBase<DerivedS>& S,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J,
Eigen::SparseMatrix<BCType>& BC) {
return igl::marching_tets(TV, TT, S, 0.0, SV, SF, J, BC);
}
// marching_tets( TV, TT, S, isovalue, SV, SF, J)
//
// performs the marching tetrahedra algorithm on a tet mesh defined by TV and
// TT with scalar values defined at each vertex in TV. The output is a
// triangle mesh approximating the isosurface coresponding to the value
// isovalue.
//
// Input:
// TV #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
// TT #tets x 4 array -- The indexes of each tet in the tetrahedral mesh
// S #tet_vertices x 1 array -- The values defined on each tet vertex
// isovalue scalar -- The isovalue of the level set we want to compute
//
// Output:
// SV #SV x 3 array -- The vertices of the output level surface mesh
// SF #SF x 3 array -- The face indexes of the output level surface mesh
// J #SF list of indices into TT revealing which tet each face comes from
template <typename DerivedTV,
typename DerivedTT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ>
IGL_INLINE void marching_tets(
const Eigen::PlainObjectBase<DerivedTV>& TV,
const Eigen::PlainObjectBase<DerivedTT>& TT,
const Eigen::PlainObjectBase<DerivedS>& S,
double isovalue,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J) {
Eigen::SparseMatrix<double> _BC;
return igl::marching_tets(TV, TT, S, isovalue, SV, SF, J, _BC);
}
// marching_tets( TV, TT, S, isovalue, SV, SF, BC)
//
// performs the marching tetrahedra algorithm on a tet mesh defined by TV and
// TT with scalar values defined at each vertex in TV. The output is a
// triangle mesh approximating the isosurface coresponding to the value
// isovalue.
//
// Input:
// TV #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
// TT #tets x 4 array -- The indexes of each tet in the tetrahedral mesh
// S #tet_vertices x 1 array -- The values defined on each tet vertex
// isovalue scalar -- The isovalue of the level set we want to compute
//
// Output:
// SV #SV x 3 array -- The vertices of the output level surface mesh
// SF #SF x 3 array -- The face indexes of the output level surface mesh
// BC #SV x #TV list of barycentric coordinates so that SV = BC*TV
template <typename DerivedTV,
typename DerivedTT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename BCType>
IGL_INLINE void marching_tets(
const Eigen::PlainObjectBase<DerivedTV>& TV,
const Eigen::PlainObjectBase<DerivedTT>& TT,
const Eigen::PlainObjectBase<DerivedS>& S,
double isovalue,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::SparseMatrix<BCType>& BC) {
Eigen::VectorXi _J;
return igl::marching_tets(TV, TT, S, isovalue, SV, SF, _J, BC);
}
// marching_tets( TV, TT, S, isovalue, SV, SF)
//
// performs the marching tetrahedra algorithm on a tet mesh defined by TV and
// TT with scalar values defined at each vertex in TV. The output is a
// triangle mesh approximating the isosurface coresponding to the value
// isovalue.
//
// Input:
// TV #tet_vertices x 3 array -- The vertices of the tetrahedral mesh
// TT #tets x 4 array -- The indexes of each tet in the tetrahedral mesh
// S #tet_vertices x 1 array -- The values defined on each tet vertex
// isovalue scalar -- The isovalue of the level set we want to compute
//
// Output:
// SV #SV x 3 array -- The vertices of the output level surface mesh
// SF #SF x 3 array -- The face indexes of the output level surface mesh
template <typename DerivedTV,
typename DerivedTT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF>
IGL_INLINE void marching_tets(
const Eigen::PlainObjectBase<DerivedTV>& TV,
const Eigen::PlainObjectBase<DerivedTT>& TT,
const Eigen::PlainObjectBase<DerivedS>& S,
double isovalue,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF) {
Eigen::VectorXi _J;
Eigen::SparseMatrix<double> _BC;
return igl::marching_tets(TV, TT, S, isovalue, SV, SF, _J, _BC);
}
}
#ifndef IGL_STATIC_LIBRARY
# include "marching_tets.cpp"
#endif
#endif // IGL_MARCHING_TETS_H
-356
View File
@@ -1,356 +0,0 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include "slice_tets.h"
#include "LinSpaced.h"
#include "sort.h"
#include "edges.h"
#include "slice.h"
#include "cat.h"
#include "ismember.h"
#include "unique_rows.h"
#include <cassert>
#include <algorithm>
#include <vector>
template <
typename DerivedV,
typename DerivedT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ,
typename BCType>
IGL_INLINE void igl::slice_tets(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedT>& T,
const Eigen::MatrixBase<DerivedS> & S,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J,
Eigen::SparseMatrix<BCType> & BC)
{
Eigen::MatrixXi sE;
Eigen::Matrix<typename DerivedSV::Scalar,Eigen::Dynamic,1> lambda;
igl::slice_tets(V,T,S,SV,SF,J,sE,lambda);
const int ns = SV.rows();
std::vector<Eigen::Triplet<BCType> > BCIJV(ns*2);
for(int i = 0;i<ns;i++)
{
BCIJV[2*i+0] = Eigen::Triplet<BCType>(i,sE(i,0), lambda(i));
BCIJV[2*i+1] = Eigen::Triplet<BCType>(i,sE(i,1),1.0-lambda(i));
}
BC.resize(SV.rows(),V.rows());
BC.setFromTriplets(BCIJV.begin(),BCIJV.end());
}
template <
typename DerivedV,
typename DerivedT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ>
IGL_INLINE void igl::slice_tets(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedT>& T,
const Eigen::MatrixBase<DerivedS> & S,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J)
{
Eigen::MatrixXi sE;
Eigen::Matrix<typename DerivedSV::Scalar,Eigen::Dynamic,1> lambda;
igl::slice_tets(V,T,S,SV,SF,J,sE,lambda);
}
template <
typename DerivedV,
typename DerivedT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ,
typename DerivedsE,
typename Derivedlambda>
IGL_INLINE void igl::slice_tets(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedT>& T,
const Eigen::MatrixBase<DerivedS> & S,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J,
Eigen::PlainObjectBase<DerivedsE>& sE,
Eigen::PlainObjectBase<Derivedlambda>& lambda)
{
using namespace Eigen;
using namespace std;
assert(V.cols() == 3 && "V should be #V by 3");
assert(T.cols() == 4 && "T should be #T by 4");
static const Eigen::Matrix<int,12,4> flipped_order =
(Eigen::Matrix<int,12,4>(12,4)<<
3,2,0,1,
3,1,2,0,
3,0,1,2,
2,3,1,0,
2,1,0,3,
2,0,3,1,
1,3,0,2,
1,2,3,0,
1,0,2,3,
0,3,2,1,
0,2,1,3,
0,1,3,2
).finished();
// number of tets
const size_t m = T.rows();
typedef typename DerivedS::Scalar Scalar;
typedef typename DerivedT::Scalar Index;
typedef Matrix<Scalar,Dynamic,1> VectorXS;
typedef Matrix<Scalar,Dynamic,4> MatrixX4S;
typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
typedef Matrix<Scalar,Dynamic,2> MatrixX2S;
typedef Matrix<Index,Dynamic,4> MatrixX4I;
typedef Matrix<Index,Dynamic,3> MatrixX3I;
typedef Matrix<Index,Dynamic,2> MatrixX2I;
typedef Matrix<Index,Dynamic,1> VectorXI;
typedef Array<bool,Dynamic,1> ArrayXb;
MatrixX4S IT(m,4);
for(size_t t = 0;t<m;t++)
{
for(size_t c = 0;c<4;c++)
{
IT(t,c) = S(T(t,c));
}
}
// Essentially, just a glorified slice(X,1)
//
// Inputs:
// T #T by 4 list of tet indices into V
// IT #IT by 4 list of isosurface values at each tet
// I #I list of bools whether to grab data corresponding to each tet
const auto & extract_rows = [](
const MatrixBase<DerivedT> & T,
const MatrixX4S & IT,
const ArrayXb & I,
MatrixX4I & TI,
MatrixX4S & ITI,
VectorXI & JI)
{
const Index num_I = std::count(I.data(),I.data()+I.size(),true);
TI.resize(num_I,4);
ITI.resize(num_I,4);
JI.resize(num_I,1);
{
size_t k = 0;
for(size_t t = 0;t<(size_t)T.rows();t++)
{
if(I(t))
{
TI.row(k) = T.row(t);
ITI.row(k) = IT.row(t);
JI(k) = t;
k++;
}
}
assert(k == num_I);
}
};
ArrayXb I13 = (IT.array()<0).rowwise().count()==1;
ArrayXb I31 = (IT.array()>0).rowwise().count()==1;
ArrayXb I22 = (IT.array()<0).rowwise().count()==2;
MatrixX4I T13,T31,T22;
MatrixX4S IT13,IT31,IT22;
VectorXI J13,J31,J22;
extract_rows(T,IT,I13,T13,IT13,J13);
extract_rows(T,IT,I31,T31,IT31,J31);
extract_rows(T,IT,I22,T22,IT22,J22);
const auto & apply_sort4 = [] (
const MatrixX4I & T,
const MatrixX4I & sJ,
MatrixX4I & sT)
{
sT.resize(T.rows(),4);
for(size_t t = 0;t<(size_t)T.rows();t++)
{
for(size_t c = 0;c<4;c++)
{
sT(t,c) = T(t,sJ(t,c));
}
}
};
const auto & apply_sort2 = [] (
const MatrixX2I & E,
const MatrixX2I & sJ,
Eigen::PlainObjectBase<DerivedsE>& sE)
{
sE.resize(E.rows(),2);
for(size_t t = 0;t<(size_t)E.rows();t++)
{
for(size_t c = 0;c<2;c++)
{
sE(t,c) = E(t,sJ(t,c));
}
}
};
const auto & one_below = [&apply_sort4](
const MatrixX4I & T,
const MatrixX4S & IT,
MatrixX2I & U,
MatrixX3I & SF)
{
// Number of tets
const size_t m = T.rows();
if(m == 0)
{
U.resize(0,2);
SF.resize(0,3);
return;
}
MatrixX4S sIT;
MatrixX4I sJ;
sort(IT,2,true,sIT,sJ);
MatrixX4I sT;
apply_sort4(T,sJ,sT);
U.resize(3*m,2);
U<<
sT.col(0),sT.col(1),
sT.col(0),sT.col(2),
sT.col(0),sT.col(3);
SF.resize(m,3);
for(size_t c = 0;c<3;c++)
{
SF.col(c) =
igl::LinSpaced<
Eigen::Matrix<typename DerivedSF::Scalar,Eigen::Dynamic,1> >
(m,0+c*m,(m-1)+c*m);
}
ArrayXb flip;
{
VectorXi _;
ismember_rows(sJ,flipped_order,flip,_);
}
for(int i = 0;i<m;i++)
{
if(flip(i))
{
SF.row(i) = SF.row(i).reverse().eval();
}
}
};
const auto & two_below = [&apply_sort4](
const MatrixX4I & T,
const MatrixX4S & IT,
MatrixX2I & U,
MatrixX3I & SF)
{
// Number of tets
const size_t m = T.rows();
if(m == 0)
{
U.resize(0,2);
SF.resize(0,3);
return;
}
MatrixX4S sIT;
MatrixX4I sJ;
sort(IT,2,true,sIT,sJ);
MatrixX4I sT;
apply_sort4(T,sJ,sT);
U.resize(4*m,2);
U<<
sT.col(0),sT.col(2),
sT.col(0),sT.col(3),
sT.col(1),sT.col(2),
sT.col(1),sT.col(3);
SF.resize(2*m,3);
SF.block(0,0,m,1) = igl::LinSpaced<VectorXI >(m,0+0*m,(m-1)+0*m);
SF.block(0,1,m,1) = igl::LinSpaced<VectorXI >(m,0+1*m,(m-1)+1*m);
SF.block(0,2,m,1) = igl::LinSpaced<VectorXI >(m,0+3*m,(m-1)+3*m);
SF.block(m,0,m,1) = igl::LinSpaced<VectorXI >(m,0+0*m,(m-1)+0*m);
SF.block(m,1,m,1) = igl::LinSpaced<VectorXI >(m,0+3*m,(m-1)+3*m);
SF.block(m,2,m,1) = igl::LinSpaced<VectorXI >(m,0+2*m,(m-1)+2*m);
ArrayXb flip;
{
VectorXi _;
ismember_rows(sJ,flipped_order,flip,_);
}
for(int i = 0;i<m;i++)
{
if(flip(i))
{
SF.row(i ) = SF.row(i ).reverse().eval();
SF.row(i+m) = SF.row(i+m).reverse().eval();
}
}
};
MatrixX3I SF13,SF31,SF22;
MatrixX2I U13,U31,U22;
one_below(T13, IT13,U13,SF13);
one_below(T31,-IT31,U31,SF31);
two_below(T22, IT22,U22,SF22);
// https://forum.kde.org/viewtopic.php?f=74&t=107974
const MatrixX2I U =
(MatrixX2I(U13.rows()+ U31.rows()+ U22.rows(),2)<<U13,U31,U22).finished();
MatrixX2I sU;
{
MatrixX2I _;
sort(U,2,true,sU,_);
}
MatrixX2I E;
VectorXI uI,uJ;
unique_rows(sU,E,uI,uJ);
MatrixX2S IE(E.rows(),2);
for(size_t t = 0;t<E.rows();t++)
{
for(size_t c = 0;c<2;c++)
{
IE(t,c) = S(E(t,c));
}
}
MatrixX2S sIE;
MatrixX2I sJ;
sort(IE,2,true,sIE,sJ);
apply_sort2(E,sJ,sE);
lambda = sIE.col(1).array() / (sIE.col(1)-sIE.col(0)).array();
SV.resize(sE.rows(),V.cols());
for(int e = 0;e<sE.rows();e++)
{
SV.row(e) = V.row(sE(e,0)).template cast<Scalar>()*lambda(e) +
V.row(sE(e,1)).template cast<Scalar>()*(1.0-lambda(e));
}
SF.resize( SF13.rows()+SF31.rows()+SF22.rows(),3);
SF<<
SF13,
U13.rows()+ SF31.rowwise().reverse().array(),
U13.rows()+U31.rows()+SF22.array();
std::for_each(
SF.data(),
SF.data()+SF.size(),
[&uJ](typename DerivedSF::Scalar & i){i=uJ(i);});
J.resize(SF.rows());
J<<J13,J31,J22,J22;
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::slice_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::Matrix<int, -1, 1, 0, -1, 1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::SparseMatrix<double, 0, int>&);
#endif
-96
View File
@@ -1,96 +0,0 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#ifndef IGL_SLICE_TETS_H
#define IGL_SLICE_TETS_H
#include "igl_inline.h"
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <vector>
namespace igl
{
// SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
// implicit equation).
//
// Inputs:
// V #V by 3 list of tet mesh vertices
// T #T by 4 list of tet indices into V
//// plane list of 4 coefficients in the plane equation: [x y z 1]'*plane = 0
// S #V list of values so that S = 0 is the desired isosurface
// Outputs:
// SV #SV by 3 list of triangle mesh vertices along slice
// SF #SF by 3 list of triangles indices into SV
// J #SF list of indices into T revealing from which tet each faces comes
// BC #SU by #V list of barycentric coordinates (or more generally: linear
// interpolation coordinates) so that SV = BC*V
//
template <
typename DerivedV,
typename DerivedT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ,
typename BCType>
IGL_INLINE void slice_tets(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedT>& T,
const Eigen::MatrixBase<DerivedS> & S,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J,
Eigen::SparseMatrix<BCType> & BC);
template <
typename DerivedV,
typename DerivedT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ>
IGL_INLINE void slice_tets(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedT>& T,
const Eigen::MatrixBase<DerivedS> & S,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J);
// Outputs:
// sE #SV by 2 list of sorted edge indices into V
// lambda #SV by 1 list of parameters along each edge in sE so that:
// SV(i,:) = V(sE(i,1),:)*lambda(i) + V(sE(i,2),:)*(1-lambda(i));
template <
typename DerivedV,
typename DerivedT,
typename DerivedS,
typename DerivedSV,
typename DerivedSF,
typename DerivedJ,
typename DerivedsE,
typename Derivedlambda
>
IGL_INLINE void slice_tets(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedT>& T,
const Eigen::MatrixBase<DerivedS> & S,
Eigen::PlainObjectBase<DerivedSV>& SV,
Eigen::PlainObjectBase<DerivedSF>& SF,
Eigen::PlainObjectBase<DerivedJ>& J,
Eigen::PlainObjectBase<DerivedsE>& sE,
Eigen::PlainObjectBase<Derivedlambda>& lambda);
}
#ifndef IGL_STATIC_LIBRARY
# include "slice_tets.cpp"
#endif
#endif
+1
View File
@@ -89,4 +89,5 @@ template bool igl::writeDMAT<Eigen::Matrix<double, 1, 3, 1, 1, 3> >(std::string,
template bool igl::writeDMAT<Eigen::Matrix<float, 1, 3, 1, 1, 3> >(std::string, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, bool);
template bool igl::writeDMAT<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool);
template bool igl::writeDMAT<Eigen::Matrix<int, -1, 2, 0, -1, 2> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, bool);
template bool igl::writeDMAT<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool);
#endif
+1 -1
View File
@@ -1275,7 +1275,7 @@ const char *__doc_igl_slice_mask = R"igl_Qu8mg5v7(// Act like the matlab X(row_m
// Y #trues-in-R by #trues-in-C matrix
//
// See also: slice_mask)igl_Qu8mg5v7";
const char *__doc_igl_slice_tets = R"igl_Qu8mg5v7(// SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
const char *__doc_igl_marching_tets = R"igl_Qu8mg5v7(// SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
// implicit equation).
//
// Inputs:
+2 -2
View File
@@ -108,7 +108,7 @@ extern const char *__doc_igl_signed_distance_winding_number;
extern const char *__doc_igl_slice;
extern const char *__doc_igl_slice_into;
extern const char *__doc_igl_slice_mask;
extern const char *__doc_igl_slice_tets;
extern const char *__doc_igl_marching_tets;
extern const char *__doc_igl_sortrows;
extern const char *__doc_igl_streamlines_init;
extern const char *__doc_igl_streamlines_next;
@@ -127,4 +127,4 @@ extern const char *__doc_igl_winding_number_2;
extern const char *__doc_igl_writeMESH;
extern const char *__doc_igl_writeOBJ;
extern const char *__doc_igl_writePLY;
extern const char *__doc_igl_readPLY;
extern const char *__doc_igl_readPLY;
+2 -2
View File
@@ -92,7 +92,7 @@
#include <igl/slice.h>
#include <igl/slice_into.h>
#include <igl/slice_mask.h>
#include <igl/slice_tets.h>
#include <igl/marching_tets.h>
#include <igl/sortrows.h>
#include <igl/triangle_triangle_adjacency.h>
#include <igl/unique.h>
@@ -191,7 +191,7 @@ void python_export_igl(py::module &m)
#include "py_igl/py_slice.cpp"
#include "py_igl/py_slice_into.cpp"
#include "py_igl/py_slice_mask.cpp"
#include "py_igl/py_slice_tets.cpp"
#include "py_igl/py_marching_tets.cpp"
#include "py_igl/py_sortrows.cpp"
#include "py_igl/py_triangle_triangle_adjacency.cpp"
#include "py_igl/py_unique.cpp"
@@ -5,7 +5,7 @@
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
m.def("slice_tets", []
m.def("marching_tets", []
(
const Eigen::MatrixXd& V,
const Eigen::MatrixXi& T,
@@ -21,8 +21,8 @@ m.def("slice_tets", []
if (plane.size() != 0)
planev = plane;
Eigen::VectorXi Jv;
igl::slice_tets(V, T, planev, U, G, Jv, BC);
igl::marching_tets(V, T, planev, U, G, Jv, BC);
J = Jv;
}, __doc_igl_slice_tets,
}, __doc_igl_marching_tets,
py::arg("V"), py::arg("T"), py::arg("plane"), py::arg("U"), py::arg("G"), py::arg("J"), py::arg("BC"));
+1 -1
View File
@@ -147,7 +147,7 @@ PYBIND11_PLUGIN(pyigl) {
slice
slice_into
slice_mask
slice_tets
marching_tets
sortrows
streamlines
triangle_triangle_adjacency
+1 -1
View File
@@ -35,7 +35,7 @@ def update(viewer):
F_vis = igl.eigen.MatrixXi()
J = igl.eigen.MatrixXi()
bary = igl.eigen.SparseMatrixd()
igl.slice_tets(V, T, plane, V_vis, F_vis, J, bary)
igl.marching_tets(V, T, plane, V_vis, F_vis, J, bary)
W_vis = igl.eigen.MatrixXd()
igl.slice(W, J, W_vis)
C_vis = igl.eigen.MatrixXd()
+1 -1
View File
@@ -55,7 +55,7 @@ def update_visualization(viewer):
# Extract triangle mesh slice through volume mesh and subdivide nasty triangles
J = igl.eigen.MatrixXi()
bary = igl.eigen.SparseMatrixd()
igl.slice_tets(V, T, plane, V_vis, F_vis, J, bary)
igl.marching_tets(V, T, plane, V_vis, F_vis, J, bary)
max_l = 0.03
while True:
l = igl.eigen.MatrixXd()
+2 -2
View File
@@ -3,7 +3,7 @@
#include <igl/parula.h>
#include <igl/readMESH.h>
#include <igl/slice.h>
#include <igl/slice_tets.h>
#include <igl/marching_tets.h>
#include <igl/winding_number.h>
#include <igl/opengl/glfw/Viewer.h>
#include <Eigen/Sparse>
@@ -40,7 +40,7 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
V.col(1)*plane(1) +
V.col(2)*plane(2)).array()
+ plane(3);
igl::slice_tets(V,T,IV,V_vis,F_vis,J,bary);
igl::marching_tets(V,T,IV,V_vis,F_vis,J,bary);
}
VectorXd W_vis;
igl::slice(W,J,W_vis);
+2 -2
View File
@@ -8,7 +8,7 @@
#include <igl/readMESH.h>
#include <igl/signed_distance.h>
#include <igl/slice_mask.h>
#include <igl/slice_tets.h>
#include <igl/marching_tets.h>
#include <igl/upsample.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/writeOBJ.h>
@@ -48,7 +48,7 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
V.col(1)*plane(1) +
V.col(2)*plane(2)).array()
+ plane(3);
igl::slice_tets(V,T,IV,V_vis,F_vis,J,bary);
igl::marching_tets(V,T,IV,V_vis,F_vis,J,bary);
igl::writeOBJ("vis.obj",V_vis,F_vis);
}
while(true)
+5
View File
@@ -0,0 +1,5 @@
get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
project(${PROJECT_NAME})
add_executable(${PROJECT_NAME}_bin main.cpp)
target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw igl::tetgen tutorials)
+48
View File
@@ -0,0 +1,48 @@
#include <igl/opengl/glfw/Viewer.h>
#include <igl/copyleft/tetgen/tetrahedralize.h>
#include <igl/readOBJ.h>
#include <igl/marching_tets.h>
#include <Eigen/Core>
#include "tutorial_shared_path.h"
int main(int argc, char * argv[])
{
// Load a surface mesh which is a cube
Eigen::MatrixXd surfaceV;
Eigen::MatrixXi surfaceF;
igl::readOBJ(TUTORIAL_SHARED_PATH "/cube.obj", surfaceV, surfaceF);
// Find the centroid of the loaded mesh
Eigen::RowVector3d surfaceCenter = surfaceV.colwise().sum() / surfaceV.rows();
// Center the mesh about the origin
surfaceV.rowwise() -= surfaceCenter;
// Tetrahedralize the surface mesh
Eigen::MatrixXd TV; // Tet mesh vertices
Eigen::MatrixXi TF; // Tet mesh boundary face indices
Eigen::MatrixXi TT; // Tet mesh tetrahedron indices
igl::copyleft::tetgen::tetrahedralize(surfaceV, surfaceF, "pq1.414a0.0001", TV, TT, TF);
// Compute a scalar at each tet vertex which is the distance from the vertex to the origin
Eigen::VectorXd S = TV.rowwise().norm();
// Compute a mesh (stored in SV, SF) representing the iso-level-set for the isovalue 0.5
Eigen::MatrixXd SV;
Eigen::MatrixXi SF;
igl::marching_tets(TV, TT, S, 0.45, SV, SF);
// Draw the mesh stored in (SV, SF)
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(SV, SF);
viewer.callback_key_down =
[&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
{
viewer.data().set_face_based(true);
return true;
};
viewer.launch();
}
+3
View File
@@ -149,6 +149,9 @@ if(TUTORIALS_CHAPTER7)
add_subdirectory("711_Subdivision")
add_subdirectory("712_DataSmoothing")
add_subdirectory("713_ShapeUp")
if(LIBIGL_WITH_TETGEN)
add_subdirectory("714_MarchingTets")
endif()
endif()