#include "outer_hull.h" #include "outer_facet.h" #include "facet_components.h" #include "winding_number.h" #include "triangle_triangle_adjacency.h" #include "unique_edge_map.h" #include "barycenter.h" #include "per_face_normals.h" #include "all_edges.h" #include "colon.h" #include "get_seconds.h" #include #include #include #include #include //#define IGL_OUTER_HULL_DEBUG template < typename DerivedV, typename DerivedF, typename DerivedG, typename DerivedJ, typename Derivedflip> IGL_INLINE void igl::outer_hull( const Eigen::PlainObjectBase & V, const Eigen::PlainObjectBase & F, Eigen::PlainObjectBase & G, Eigen::PlainObjectBase & J, Eigen::PlainObjectBase & flip) { using namespace Eigen; using namespace std; using namespace igl; typedef typename DerivedF::Index Index; Matrix C; typedef Matrix MatrixXV; typedef Matrix MatrixXF; typedef Matrix MatrixXG; typedef Matrix MatrixXJ; const Index m = F.rows(); #ifdef IGL_OUTER_HULL_DEBUG cout<<"outer hull..."< MatrixX2I; typedef Matrix VectorXI; MatrixX2I E,uE; VectorXI EMAP; vector > uE2E; unique_edge_map(F,E,uE,EMAP,uE2E); vector > > TT,_1; triangle_triangle_adjacency(E,EMAP,uE2E,false,TT,_1); VectorXI counts; #ifdef IGL_OUTER_HULL_DEBUG cout<<"facet components..."< N; #ifdef IGL_OUTER_HULL_DEBUG cout<<"normals..."< FH(m,false); vector EH(3*m,false); vector vG(ncc); vector vJ(ncc); vector vIM(ncc); for(size_t id = 0;id g(ncc,0); // place order of each face in its respective component for(Index f = 0;f Q; Q.push(f+0*m); Q.push(f+1*m); Q.push(f+2*m); flip(f) = f_flip; #ifdef IGL_OUTER_HULL_DEBUG cout<<"BFS..."< max_nN; for(const auto & ne : neighbors) { const int nf = ne%m; if(nf == f) { continue; } const int nc = ne/m; // are faces consistently oriented //const int ns = F(nf,(nc+1)%3); const int nd = F(nf,(nc+2)%3); const bool cons = (flip(f)?fd:fs) == nd; const auto & nN = (cons? (flip(f)?-1:1.) : (flip(f)?1.:-1.) )*N.row(nf); const auto & ndi = M_PI - atan2( fN.cross(nN).dot(eV), fN.dot(nN)); if(ndi>=max_di) { max_ne = ne; max_di = ndi; max_nN = nN; } } if(max_ne>=0) { const int nf = max_ne%m; const int nc = max_ne/m; const int nd = F(nf,(nc+2)%3); const bool cons = (flip(f)?fd:fs) == nd; flip(nf) = (cons ? flip(f) : !flip(f)); const int ne1 = nf+((nc+1)%3)*m; const int ne2 = nf+((nc+2)%3)*m; if(!EH[ne1]) { Q.push(ne1); } if(!EH[ne2]) { Q.push(ne2); } } } { vG[id].resize(FHcount,3); vJ[id].resize(FHcount,1); //nG += FHcount; size_t h = 0; assert(counts(id) == IM.rows()); for(int i = 0;i & V, const MatrixXV & BC, const MatrixXG & A, const MatrixXJ & AJ, const MatrixXG & B)->bool { const auto & bounding_box = []( const Eigen::PlainObjectBase & V, const MatrixXG & F)-> MatrixXV { MatrixXV BB(2,3); BB<< 1e26,1e26,1e26, -1e26,-1e26,-1e26; const size_t m = F.rows(); for(size_t f = 0;f0 || (ABB.row(0)-BBB.row(1)).maxCoeff()>0 ) { // bounding boxes do not overlap return false; } //////////////////////////////////////////////////////////////////////// // POTENTIAL ROBUSTNESS WEAK AREA //////////////////////////////////////////////////////////////////////// // // q could be so close (<~1e-16) to B that the winding number is not a robust way to // determine inside/outsideness. We could try to find a _better_ q which is // farther away, but couldn't they all be bad? MatrixXV q = BC.row(AJ(0)); // In a perfect world, it's enough to test a single point. double w; winding_number_3( V.data(),V.rows(), B.data(),B.rows(), q.data(),1,&w); return fabs(w)>0.5; }; // Reject components which are completely inside other components vector keep(ncc,true); size_t nG = 0; // This is O( ncc * ncc * m) for(size_t id = 0;id, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); #endif