Adding an AABB tree to EigenMesh3D.

Yet to be used.
This commit is contained in:
tamasmeszaros 2019-01-14 17:28:02 +01:00
parent 6ac54896fa
commit 7fa430c56d
6 changed files with 155 additions and 137 deletions

View file

@ -19,6 +19,10 @@
namespace Slic3r {
namespace sla {
/* **************************************************************************
* SpatIndex implementation
* ************************************************************************** */
class SpatIndex::Impl {
public:
using BoostIndex = boost::geometry::index::rtree< SpatElement,
@ -78,6 +82,73 @@ size_t SpatIndex::size() const
return m_impl->m_store.size();
}
/* ****************************************************************************
* EigenMesh3D implementation
* ****************************************************************************/
class EigenMesh3D::AABBImpl: public igl::AABB<Eigen::MatrixXd, 3> {};
EigenMesh3D::EigenMesh3D(): m_aabb(new AABBImpl()) {}
EigenMesh3D::EigenMesh3D(const TriangleMesh& tmesh): m_aabb(new AABBImpl()) {
static const double dEPS = 1e-6;
const stl_file& stl = tmesh.stl;
auto&& bb = tmesh.bounding_box();
m_ground_level += bb.min(Z);
Eigen::MatrixXd V;
Eigen::MatrixXi F;
V.resize(3*stl.stats.number_of_facets, 3);
F.resize(stl.stats.number_of_facets, 3);
for (unsigned int i = 0; i < stl.stats.number_of_facets; ++i) {
const stl_facet* facet = stl.facet_start+i;
V(3*i+0, 0) = double(facet->vertex[0](0));
V(3*i+0, 1) = double(facet->vertex[0](1));
V(3*i+0, 2) = double(facet->vertex[0](2));
V(3*i+1, 0) = double(facet->vertex[1](0));
V(3*i+1, 1) = double(facet->vertex[1](1));
V(3*i+1, 2) = double(facet->vertex[1](2));
V(3*i+2, 0) = double(facet->vertex[2](0));
V(3*i+2, 1) = double(facet->vertex[2](1));
V(3*i+2, 2) = double(facet->vertex[2](2));
F(i, 0) = int(3*i+0);
F(i, 1) = int(3*i+1);
F(i, 2) = int(3*i+2);
}
// We will convert this to a proper 3d mesh with no duplicate points.
Eigen::VectorXi SVI, SVJ;
igl::remove_duplicate_vertices(V, F, dEPS, m_V, SVI, SVJ, m_F);
// Build the AABB accelaration tree
m_aabb->init(m_V, m_F);
}
EigenMesh3D::~EigenMesh3D() {}
EigenMesh3D::EigenMesh3D(const EigenMesh3D &other):
m_V(other.m_V), m_F(other.m_F), m_ground_level(other.m_ground_level),
m_aabb( new AABBImpl(*other.m_aabb) ) {}
EigenMesh3D &EigenMesh3D::operator=(const EigenMesh3D &other)
{
m_V = other.m_V;
m_F = other.m_F;
m_ground_level = other.m_ground_level;
m_aabb.reset(new AABBImpl(*other.m_aabb)); return *this;
}
double EigenMesh3D::query_ray_hit(const Vec3d &s, const Vec3d &dir) const
{
return 0;
}
bool point_on_edge(const Vec3d& p, const Vec3d& e1, const Vec3d& e2,
double eps = 0.05)
{
@ -93,10 +164,10 @@ template<class Vec> double distance(const Vec& pp1, const Vec& pp2) {
return std::sqrt(p.transpose() * p);
}
PointSet normals(const PointSet& points, const EigenMesh3D& emesh,
PointSet normals(const PointSet& points, const EigenMesh3D& mesh,
double eps,
std::function<void()> throw_on_cancel) {
if(points.rows() == 0 || emesh.V.rows() == 0 || emesh.F.rows() == 0)
if(points.rows() == 0 || mesh.V().rows() == 0 || mesh.F().rows() == 0)
return {};
Eigen::VectorXd dists;
@ -105,23 +176,24 @@ PointSet normals(const PointSet& points, const EigenMesh3D& emesh,
// We need to remove duplicate vertices and have a true index triangle
// structure
/*
EigenMesh3D mesh;
Eigen::VectorXi SVI, SVJ;
static const double dEPS = 1e-6;
igl::remove_duplicate_vertices(emesh.V, emesh.F, dEPS,
mesh.V, SVI, SVJ, mesh.F);
mesh.V, SVI, SVJ, mesh.F);*/
igl::point_mesh_squared_distance( points, mesh.V, mesh.F, dists, I, C);
igl::point_mesh_squared_distance( points, mesh.V(), mesh.F(), dists, I, C);
PointSet ret(I.rows(), 3);
for(int i = 0; i < I.rows(); i++) {
throw_on_cancel();
auto idx = I(i);
auto trindex = mesh.F.row(idx);
auto trindex = mesh.F().row(idx);
const Vec3d& p1 = mesh.V.row(trindex(0));
const Vec3d& p2 = mesh.V.row(trindex(1));
const Vec3d& p3 = mesh.V.row(trindex(2));
const Vec3d& p1 = mesh.V().row(trindex(0));
const Vec3d& p2 = mesh.V().row(trindex(1));
const Vec3d& p3 = mesh.V().row(trindex(2));
// We should check if the point lies on an edge of the hosting triangle.
// If it does than all the other triangles using the same two points
@ -159,18 +231,18 @@ PointSet normals(const PointSet& points, const EigenMesh3D& emesh,
// vector for the neigboring triangles including the detected one.
std::vector<Vec3i> neigh;
if(ic >= 0) { // The point is right on a vertex of the triangle
for(int n = 0; n < mesh.F.rows(); ++n) {
for(int n = 0; n < mesh.F().rows(); ++n) {
throw_on_cancel();
Vec3i ni = mesh.F.row(n);
Vec3i ni = mesh.F().row(n);
if((ni(X) == ic || ni(Y) == ic || ni(Z) == ic))
neigh.emplace_back(ni);
}
}
else if(ia >= 0 && ib >= 0) { // the point is on and edge
// now get all the neigboring triangles
for(int n = 0; n < mesh.F.rows(); ++n) {
for(int n = 0; n < mesh.F().rows(); ++n) {
throw_on_cancel();
Vec3i ni = mesh.F.row(n);
Vec3i ni = mesh.F().row(n);
if((ni(X) == ia || ni(Y) == ia || ni(Z) == ia) &&
(ni(X) == ib || ni(Y) == ib || ni(Z) == ib))
neigh.emplace_back(ni);
@ -180,9 +252,9 @@ PointSet normals(const PointSet& points, const EigenMesh3D& emesh,
// Calculate the normals for the neighboring triangles
std::vector<Vec3d> neighnorms; neighnorms.reserve(neigh.size());
for(const Vec3i& tri : neigh) {
const Vec3d& pt1 = mesh.V.row(tri(0));
const Vec3d& pt2 = mesh.V.row(tri(1));
const Vec3d& pt3 = mesh.V.row(tri(2));
const Vec3d& pt1 = mesh.V().row(tri(0));
const Vec3d& pt2 = mesh.V().row(tri(1));
const Vec3d& pt3 = mesh.V().row(tri(2));
Eigen::Vector3d U = pt2 - pt1;
Eigen::Vector3d V = pt3 - pt1;
neighnorms.emplace_back(U.cross(V).normalized());
@ -228,10 +300,24 @@ double ray_mesh_intersect(const Vec3d& s,
{
igl::Hit hit;
hit.t = std::numeric_limits<float>::infinity();
igl::ray_mesh_intersect(s, dir, m.V, m.F, hit);
// Fck: this does not use any kind of spatial index acceleration...
igl::ray_mesh_intersect(s, dir, m.V(), m.F(), hit);
return double(hit.t);
}
// An enhanced version of ray_mesh_intersect for the pinheads. This will shoot
// multiple rays to detect collisions more accurately.
double pinhead_mesh_intersect(const Vec3d& jp,
const Vec3d& dir,
double r1,
double r2,
const EigenMesh3D& m)
{
return 0;
}
// Clustering a set of points by the given criteria
ClusteredPoints cluster(
const sla::PointSet& points,
@ -309,53 +395,5 @@ ClusteredPoints cluster(
return result;
}
using Segments = std::vector<std::pair<Vec2d, Vec2d>>;
Segments model_boundary(const EigenMesh3D& emesh, double offs)
{
Segments ret;
Polygons pp;
pp.reserve(size_t(emesh.F.rows()));
for (int i = 0; i < emesh.F.rows(); i++) {
auto trindex = emesh.F.row(i);
auto& p1 = emesh.V.row(trindex(0));
auto& p2 = emesh.V.row(trindex(1));
auto& p3 = emesh.V.row(trindex(2));
Polygon p;
p.points.resize(3);
p.points[0] = Point::new_scale(p1(X), p1(Y));
p.points[1] = Point::new_scale(p2(X), p2(Y));
p.points[2] = Point::new_scale(p3(X), p3(Y));
p.make_counter_clockwise();
pp.emplace_back(p);
}
ExPolygons merged = union_ex(Slic3r::offset(pp, float(scale_(offs))), true);
for(auto& expoly : merged) {
auto lines = expoly.lines();
for(Line& l : lines) {
Vec2d a(l.a(X) * SCALING_FACTOR, l.a(Y) * SCALING_FACTOR);
Vec2d b(l.b(X) * SCALING_FACTOR, l.b(Y) * SCALING_FACTOR);
ret.emplace_back(std::make_pair(a, b));
}
}
return ret;
}
//struct SegmentIndex {
//};
//using SegmentIndexEl = std::pair<Segment, unsigned>;
//SegmentIndexEl
}
}