// This file is part of libigl, a simple c++ geometry processing library. // // Copyright (C) 2015 Alec Jacobson // // 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 "remesh_intersections.h" #include "SelfIntersectMesh.h" #include "assign_scalar.h" #include "projected_delaunay.h" #include #include template < typename DerivedV, typename DerivedF, typename Kernel, typename DerivedVV, typename DerivedFF, typename DerivedJ, typename DerivedIM> IGL_INLINE void igl::cgal::remesh_intersections( const Eigen::PlainObjectBase & V, const Eigen::PlainObjectBase & F, const std::vector > & T, const std::map< typename DerivedF::Index, std::pair > > & offending, const std::map< std::pair, std::vector > & edge2faces, Eigen::PlainObjectBase & VV, Eigen::PlainObjectBase & FF, Eigen::PlainObjectBase & J, Eigen::PlainObjectBase & IM) { using namespace std; using namespace Eigen; typedef typename DerivedF::Index Index; typedef CGAL::Point_3 Point_3; //typedef CGAL::Segment_3 Segment_3; //typedef CGAL::Triangle_3 Triangle_3; typedef CGAL::Plane_3 Plane_3; //typedef CGAL::Point_2 Point_2; //typedef CGAL::Segment_2 Segment_2; //typedef CGAL::Triangle_2 Triangle_2; typedef CGAL::Triangulation_vertex_base_2 TVB_2; typedef CGAL::Constrained_triangulation_face_base_2 CTFB_2; typedef CGAL::Triangulation_data_structure_2 TDS_2; typedef CGAL::Exact_intersections_tag Itag; typedef CGAL::Constrained_Delaunay_triangulation_2 CDT_2; typedef CGAL::Constrained_triangulation_plus_2 CDT_plus_2; typedef std::pair EMK; typedef std::vector EMV; //typedef std::map EdgeMap; typedef std::pair EMK; //typedef std::vector ObjectList; typedef std::vector IndexList; int NF_count = 0; // list of new faces, we'll fix F later vector< typename Eigen::Matrix > NF(offending.size()); // list of new vertices typedef vector Point_3List; Point_3List NV; Index NV_count = 0; vector cdt(offending.size()); vector P(offending.size()); // Use map for *all* faces map v2i; #ifdef IGL_SELFINTERSECTMESH_DEBUG double t_proj_del = 0; #endif // Unfortunately it looks like CGAL has trouble allocating memory when // multiple openmp threads are running. Crashes durring CDT... //// Loop over offending triangles //const size_t noff = offending.size(); //# pragma omp parallel for if (noff>1000) for(const auto & okv : offending) { // index in F const Index f = okv.first; const Index o = okv.second.first; { #ifdef IGL_SELFINTERSECTMESH_DEBUG const double t_before = get_seconds(); #endif CDT_plus_2 cdt_o; projected_delaunay(T[f],okv.second.second,cdt_o); cdt[o] = cdt_o; #ifdef IGL_SELFINTERSECTMESH_DEBUG t_proj_del += (get_seconds()-t_before); #endif } // Q: Is this also delaunay in 3D? // A: No, because the projection is affine and delaunay is not affine // invariant. // Q: Then, can't we first get the 2D delaunay triangulation, then lift it // to 3D and flip any offending edges? // Plane of projection (also used by projected_delaunay) P[o] = Plane_3(T[f].vertex(0),T[f].vertex(1),T[f].vertex(2)); // Build index map { Index i=0; for( typename CDT_plus_2::Finite_vertices_iterator vit = cdt[o].finite_vertices_begin(); vit != cdt[o].finite_vertices_end(); ++vit) { if(i<3) { //cout<point())?" == ":" != ")<< // P[o].to_3d(vit->point())<point())); #endif // For first three, use original index in F //# pragma omp critical v2i[vit] = F(f,i); }else { const Point_3 vit_point_3 = P[o].to_3d(vit->point()); // First look up each edge's neighbors to see if exact point has // already been added (This makes everything a bit quadratic) bool found = false; for(int e = 0; e<3 && !found;e++) { // Index of F's eth edge in V Index i = F(f,(e+1)%3); Index j = F(f,(e+2)%3); // Be sure that ij) { swap(i,j); } assert(edge2faces.count(EMK(i,j))==1); const EMV & facesij = edge2faces.find(EMK(i,j))->second; // loop over neighbors for( typename IndexList::const_iterator nit = facesij.begin(); nit != facesij.end() && !found; nit++) { // don't consider self if(*nit == f) { continue; } // index of neighbor in offending (to find its cdt) assert(offending.count(*nit) == 1); Index no = offending.find(*nit)->second.first; // Loop over vertices of that neighbor's cdt (might not have been // processed yet, but then it's OK because it'll just be empty) for( typename CDT_plus_2::Finite_vertices_iterator uit = cdt[no].finite_vertices_begin(); uit != cdt[no].finite_vertices_end() && !found; ++uit) { if(vit_point_3 == P[no].to_3d(uit->point())) { assert(v2i.count(uit) == 1); //# pragma omp critical v2i[vit] = v2i[uit]; found = true; } } } } if(!found) { //# pragma omp critical { v2i[vit] = V.rows()+NV_count; NV.push_back(vit_point_3); NV_count++; } } } i++; } } { Index i = 0; // Resize to fit new number of triangles NF[o].resize(cdt[o].number_of_faces(),3); //# pragma omp atomic NF_count+=NF[o].rows(); // Append new faces to NF for( typename CDT_plus_2::Finite_faces_iterator fit = cdt[o].finite_faces_begin(); fit != cdt[o].finite_faces_end(); ++fit) { NF[o](i,0) = v2i[fit->vertex(0)]; NF[o](i,1) = v2i[fit->vertex(1)]; NF[o](i,2) = v2i[fit->vertex(2)]; i++; } } } #ifdef IGL_SELFINTERSECTMESH_DEBUG cout<<"CDT: "<(); { Index i = 0; for( typename Point_3List::const_iterator nvit = NV.begin(); nvit != NV.end(); nvit++) { for(Index d = 0;d<3;d++) { const Point_3 & p = *nvit; // Don't convert via double if output type is same as Kernel assign_scalar(p[d], VV(V.rows()+i,d)); } i++; } } IM.resize(VV.rows(),1); map vv2i; // Safe to check for duplicates using double for original vertices: if // incoming reps are different then the points are unique. for(Index v = 0;v, Eigen::Matrix, CGAL::Epeck, Eigen::Matrix, -1, 3, 0, -1, 3>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epeck, Eigen::Matrix, -1, 3, 0, -1, 3>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epeck, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epeck, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epick, Eigen::Matrix, -1, 3, 0, -1, 3>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epick, Eigen::Matrix, -1, 3, 0, -1, 3>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epick, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epick, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epeck, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epick, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, -1, 3, 0, -1, 3>, Eigen::Matrix, CGAL::Epeck, Eigen::Matrix, -1, 3, 0, -1, 3>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, -1, 3, 0, -1, 3>, Eigen::Matrix, CGAL::Epick, Eigen::Matrix, -1, 3, 0, -1, 3>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epeck, Eigen::Matrix, -1, 3, 0, -1, 3>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::cgal::remesh_intersections, Eigen::Matrix, CGAL::Epick, Eigen::Matrix, -1, 3, 0, -1, 3>, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, std::vector, std::allocator > > const&, std::map::Index, std::pair::Index, std::vector > >, std::less::Index>, std::allocator::Index const, std::pair::Index, std::vector > > > > > const&, std::map::Index, Eigen::Matrix::Index>, std::vector::Index, std::allocator::Index> >, std::less::Index, Eigen::Matrix::Index> >, std::allocator::Index, Eigen::Matrix::Index> const, std::vector::Index, std::allocator::Index> > > > > const&, Eigen::PlainObjectBase, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); #endif