XCore intersectMesh: better closed edges extraction. OK with EPFL case.

Cassiopee/XCore/XCore/intersectMesh/dcel.cpp

@@ -358,6 +358,7 @@ void Dcel::init_hedges_and_faces(Smesh &M, E_Int color)
         Vertex *P = xit->key;
 
         Hedge *h = new Hedge(P);
+        h->eid = i;
 
         list[p].push_back(h);
 
@@ -367,6 +368,7 @@ void Dcel::init_hedges_and_faces(Smesh &M, E_Int color)
         Vertex *V = xit->key;
 
         Hedge *t = new Hedge(V);
+        t->eid = i;
 
         list[q].push_back(t);
 
@@ -392,7 +394,7 @@ void Dcel::init_hedges_and_faces(Smesh &M, E_Int color)
         const E_Float *N = &pnormals[3*pid];
         assert(Sign(K_MATH::norm(N, 3)-1) == 0);
 
-        sort_leaving_hedges(hedges, N);
+        sort_leaving_hedges(hedges, N, M);
 
         for (size_t i = 0; i < hedges.size(); i++) {
             Hedge *h = hedges[i];
@@ -419,23 +421,10 @@ void Dcel::init_hedges_and_faces(Smesh &M, E_Int color)
 
         Face *f = new Face;
         f->oid[color] = i;
-        //f->oid[color] = M.l2gf.at(i);
-
 
         assert(M.E2F[first_edge][0] == (E_Int)i || M.E2F[first_edge][1] == E_Int(i));
         Hedge *REP = (M.E2F[first_edge][0] == (E_Int)i) ? h : t;
 
-        /*
-        if (REP->left != NULL) {
-            Vertex *v = REP->orig;
-            point_write("point", v->x, v->y, v->z);
-            hedge_write("hedge", REP);
-            hedge_write("twin", REP->twin);
-            hedge_write("next", REP->next);
-            hedge_write("prev", REP->prev);
-        }
-        */
-
         assert(REP->left == NULL);
 
         f->rep = REP;
@@ -1141,7 +1130,9 @@ void Dcel::find_intersections_3D(const Smesh &M, const Smesh &S)
     }
 }
 
-void Dcel::sort_leaving_hedges(std::vector<Hedge *> &leaving, const E_Float N[3]) const
+void Dcel::sort_leaving_hedges(std::vector<Hedge *> &leaving,
+    const E_Float N[3],
+    const Smesh &M) const
 {
     // Choose a vector that is not parallel to N
 
@@ -1221,7 +1212,9 @@ void Dcel::sort_leaving_hedges(std::vector<Hedge *> &leaving, const E_Float N[3]
         else {
             Hedge *h = leaving[i];
             Hedge *w = leaving[j];
+
             assert(h->color != w->color);
+
             Vertex *P = h->orig;
             Vertex *Q = h->twin->orig;
             if (cmp_vtx(P, Q) < 0) return true;
@@ -1370,7 +1363,7 @@ void Dcel::resolve_hedges(const Smesh &M, const Smesh &S)
         E_Float NORM = K_MATH::norm(N, 3);
         assert(Sign(NORM -1) == 0);
 
-        sort_leaving_hedges(leaving, N);
+        sort_leaving_hedges(leaving, N, M);
 
         for (size_t i = 0; i < leaving.size(); i++) {
             Hedge *h = leaving[i];

Cassiopee/XCore/XCore/intersectMesh/dcel.h

@@ -108,5 +108,6 @@ struct Dcel {
 
     void trace_hedge(Hedge *sh, const Smesh &M, const Smesh &S, E_Int hid);
 
-    void sort_leaving_hedges(std::vector<Hedge *> &leaving, const E_Float N[3]) const;
+    void sort_leaving_hedges(std::vector<Hedge *> &leaving, const E_Float N[3],
+        const Smesh &M) const;
 };

Cassiopee/XCore/XCore/intersectMesh/hedge.h

@@ -36,6 +36,8 @@ struct Hedge {
     E_Int color;
     Cycle *cycle;
 
+    E_Int eid;
+
     // Projection of orig
     E_Float proj_ox = -10000;
     E_Float proj_oy = -10000;

Cassiopee/XCore/XCore/intersectMesh/intersectMesh.cpp

@@ -314,6 +314,11 @@ PyObject *K_XCORE::intersectMesh(PyObject *self, PyObject *args)
     IMesh S(*sarray.cn, sarray.X, sarray.Y, sarray.Z, sarray.npts);
 
     M.make_skin();
+    S.make_skin();
+
+    M.orient_skin(OUT);
+    S.orient_skin(IN);
+
     assert(M.patch.empty());
     for (E_Int fid : M.skin) M.patch.insert(fid);
 
@@ -334,15 +339,12 @@ PyObject *K_XCORE::intersectMesh(PyObject *self, PyObject *args)
 
     for (E_Int i = 0; i < spatch_size; i++) S.patch.insert(spatch[i]-1);
 
-    M.orient_skin(OUT);
-    S.orient_skin(IN);
-
     // Extract surface meshes
     Smesh Mf(M);
     Smesh Sf(S);
 
-    //Mf.write_ngon("Mf");
-    //Sf.write_ngon("Sf");
+    Mf.write_ngon("Mf");
+    Sf.write_ngon("Sf");
 
     puts("Making point edges...");
     Mf.make_point_edges();

Cassiopee/XCore/XCore/intersectMesh/mesh.cpp

@@ -856,6 +856,37 @@ E_Int IMesh::face_contains_point(E_Int face, E_Float x, E_Float y, E_Float z) co
     return hit;
 }
 
+void IMesh::extract_edge_points(E_Int a, E_Int b, std::list<E_Int> &points)
+{
+    E_Int ref = 0;
+
+    points.clear();
+    points.push_back(a);
+    points.push_back(b);
+
+    do {
+        ref = 0;
+
+        auto pos = points.begin();
+
+        assert(*std::prev(points.end()) == b);
+
+        for (auto it = points.begin(); it != std::prev(points.end()); it++) {
+            E_Int a = *it;
+            E_Int b = *std::next(it);
+
+            UEdge e(a, b);
+
+            auto search = ecenter.find(e);
+
+            if (search != ecenter.end()) {
+                points.insert(std::next(it), search->second);
+                ref = 1;
+            }
+        }
+    } while (ref);
+}
+
 IMesh IMesh::extract_conformized()
 {
     // Keep all the points
@@ -880,6 +911,16 @@ IMesh IMesh::extract_conformized()
             E_Int p = pn[j];
             E_Int q = pn[(j+1)%pn.size()];
 
+            std::list<E_Int> epoints;
+
+            extract_edge_points(p, q, epoints);
+
+            epoints.pop_back();
+
+            for (auto it = epoints.begin(); it != epoints.end(); it++)
+                new_face.push_back(*it);
+
+            /*
             UEdge e(p, q);
 
             auto it = ecenter.find(e);
@@ -890,6 +931,7 @@ IMesh IMesh::extract_conformized()
             } else {
                 new_face.push_back(p);
             }
+            */
         }
 
         new_nf++;

Cassiopee/XCore/XCore/intersectMesh/mesh.h

@@ -22,6 +22,7 @@
 #include <array>
 #include <map>
 #include <set>
+#include <list>
 
 #include "point.h"
 #include "xcore.h"
@@ -112,6 +113,8 @@ struct IMesh {
 
     void make_edges();
 
+    void extract_edge_points(E_Int a, E_Int b, std::list<E_Int> &pts);
+
     inline bool face_is_quad(E_Int face) const { return F[face].size() == 4; }
 
     inline bool face_is_tri(E_Int face) const { return F[face].size() == 3; }

Cassiopee/XCore/XCore/intersectMesh/prepareMeshesForIntersection.cpp

@@ -56,7 +56,7 @@ PyObject *K_XCORE::prepareMeshesForIntersection(PyObject *self, PyObject *args)
 
     puts("Preparing meshes for intersection...");
 
-    // Init and orient master/slave meshes
+    // Init master/slave meshes
     IMesh M(*marray.cn, marray.X, marray.Y, marray.Z, marray.npts);
     IMesh S(*sarray.cn, sarray.X, sarray.Y, sarray.Z, sarray.npts);
 
@@ -83,11 +83,6 @@ PyObject *K_XCORE::prepareMeshesForIntersection(PyObject *self, PyObject *args)
     //for (E_Int i = 0; i < S.nf; i++) {
         const auto &pn = S.F[fid];
         size_t stride = pn.size();
-        if (stride > 4) {
-            std::vector<Point> points;
-            for (auto p : pn) points.push_back(Point(S.X[p], S.Y[p], S.Z[p]));
-            point_write("face", points);
-        }
         assert(stride == 3 || stride == 4);
 
         E_Int keep = 1;

Cassiopee/XCore/XCore/intersectMesh/smesh.cpp

@@ -20,6 +20,7 @@
 #include "primitives.h"
 #include "triangle.h"
 #include "mesh.h"
+#include "io.h"
 
 #include <cassert>
 #include <cstdio>
@@ -51,6 +52,7 @@ bool Smesh::ccw_oriented(E_Int face)
 Smesh::Smesh(const IMesh &M)
 {
     F.resize(M.patch.size());
+    //assert(M.patch.size() == M.skin.size());
 
     nf = 0;
     np = 0;
@@ -78,8 +80,11 @@ Smesh::Smesh(const IMesh &M)
         nf++;
     }
 
+    assert(np == g2lp.size());
+    assert(np == l2gp.size());
+
     // Get the points
-    np = g2lp.size();
+    //np = g2lp.size();
     X.resize(np);
     Y.resize(np);
     Z.resize(np);
@@ -115,41 +120,51 @@ void Smesh::make_edges()
     F2E.resize(F.size());
     std::map<o_edge, E_Int, o_edge_cmp> edges;
 
+    ne = 0;
+
     for (E_Int i = 0; i < nf; i++) {
         auto &face = F[i];
         for (size_t j = 0; j < face.size(); j++) {
             E_Int p = face[j];
             E_Int q = face[(j+1)%face.size()];
-            o_edge e(p, q);
-            auto it = edges.find(e);
+            o_edge EDGE(p, q);
+            auto it = edges.find(EDGE);
             if (it == edges.end()) {
-                F2E[i].push_back(E.size());
-                edges[e] = E.size();
-                E.push_back(e);
+                F2E[i].push_back(ne);
+                edges[EDGE] = ne;
+                E.push_back(EDGE);
+                ne++;
             } else {
                 F2E[i].push_back(it->second);
             }
         }
     }
 
-    ne = E.size();
-
-    //assert(F.size()+1 + X.size() == E.size() + 2);
+    assert(ne == E.size());
 
     // Make edge_to_face
     E2F.resize(ne, {-1,-1});
+
     std::vector<E_Int> count(ne, 0);
+
     for (E_Int i = 0; i < nf; i++) {
-        auto &face = F2E[i];
+        const auto &face = F2E[i];
+        const auto &pn = F[i];
+        assert(face.size() == pn.size());
+
         for (size_t j = 0; j < face.size(); j++) {
             E_Int e = face[j];
             count[e]++;
-            assert(count[e] <= 2);
+
             if (E2F[e][0] == -1) E2F[e][0] = i;
-            else E2F[e][1] = i;
+            else {
+                assert(E2F[e][1] == -1);
+                E2F[e][1] = i;
+            }
         }
     }
 
+
     // Check
     for (E_Int i = 0; i < ne; i++) {
         E_Int fi = E2F[i][0];
@@ -739,7 +754,7 @@ Smesh Smesh::extract_conformized()
     return ret;
 }
 
-void Smesh::write_faces(const char *fname, const std::vector<E_Int> &faces)
+void Smesh::write_faces(const char *fname, const std::vector<E_Int> &faces) const
 {
     std::map<E_Int, E_Int> pmap;
 

Cassiopee/XCore/XCore/intersectMesh/smesh.h

@@ -120,7 +120,7 @@ struct Smesh {
 
     //std::vector<pointFace> locate(E_Float x, E_Float y, E_Float z) const;
 
-    void write_faces(const char *fname, const std::vector<E_Int> &faces);
+    void write_faces(const char *fname, const std::vector<E_Int> &faces) const;
 
     void write_ngon(const char *fname);