Fixed a bug in determining if outer facet is flipped.

Former-commit-id: d3df3e1db5d96b65c4833949d90d84ecacf775d7

include/igl/cgal/order_facets_around_edge.cpp

@@ -81,8 +81,22 @@ void igl::cgal::order_facets_around_edge(
                 order(1, 0) = 1;
                 break;
             case CGAL::COPLANAR:
-                order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
-                order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
+                {
+                    Plane_3 P1(ps, pd, p1);
+                    Plane_3 P2(ps, pd, p2);
+                    if (P1.orthogonal_direction() == P2.orthogonal_direction()){
+                        // Duplicated face, use index to break tie.
+                        order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
+                        order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
+                    } else {
+                        // Coplanar faces, one on each side of the edge.
+                        // It is equally valid to order them (0, 1) or (1, 0).
+                        // I cannot think of any reason to prefer one to the
+                        // other.  So just use (0, 1) ordering by default.
+                        order(0, 0) = 0;
+                        order(1, 0) = 1;
+                    }
+                }
                 break;
             default:
                 assert(false);
@@ -199,7 +213,18 @@ void igl::cgal::order_facets_around_edge(
         const Point_3& p_a = opposite_vertices[order(i, 0)];
         const Point_3& p_b =
             opposite_vertices[order((i+1)%num_adj_faces, 0)];
-        if (CGAL::orientation(p_s, p_d, p_a, p_b) == CGAL::POSITIVE) {
+        auto orientation = CGAL::orientation(p_s, p_d, p_a, p_b);
+        if (orientation == CGAL::POSITIVE) {
+            // Angle between triangle (p_s, p_d, p_a) and (p_s, p_d, p_b) is
+            // more than 180 degrees.
+            start_idx = (i+1)%num_adj_faces;
+            break;
+        } else if (orientation == CGAL::COPLANAR &&
+                Plane_3(p_s, p_d, p_a).orthogonal_direction() !=
+                Plane_3(p_s, p_d, p_b).orthogonal_direction()) {
+            // All 4 points are coplanar, but p_a and p_b are on each side of
+            // the edge (p_s, p_d).  This means the angle between triangle
+            // (p_s, p_d, p_a) and (p_s, p_d, p_b) is exactly 180 degrees.
             start_idx = (i+1)%num_adj_faces;
             break;
         }

include/igl/cgal/outer_facet.cpp

@@ -10,6 +10,7 @@
 #include "../outer_element.h"
 #include "order_facets_around_edge.h"
 #include <algorithm>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
 template<
     typename DerivedV,
@@ -36,10 +37,13 @@ IGL_INLINE void igl::cgal::outer_facet(
     //
     //    3. Because the vertex s is an outer vertex by construction (see
     //    implemnetation of outer_edge()). The first adjacent facet is facing
-    //    outside (i.e. flipped=false) if it contains directed edge (s, d).  
-    //
-    typedef typename DerivedV::Scalar Scalar; typedef typename DerivedV::Index
-        Index; const size_t INVALID = std::numeric_limits<size_t>::max();
+    //    outside (i.e. flipped=false) if it has positive X normal component.
+    //    If it has zero normal component, it is facing outside if it contains
+    //    directed edge (s, d).  
+
+    typedef typename DerivedV::Scalar Scalar;
+    typedef typename DerivedV::Index Index;
+    const size_t INVALID = std::numeric_limits<size_t>::max();
 
     Index s,d;
     Eigen::Matrix<Index,Eigen::Dynamic,1> incident_faces;
@@ -70,5 +74,18 @@ IGL_INLINE void igl::cgal::outer_facet(
     order_facets_around_edge(V, F, s, d, adj_faces, order);
 
     f = signed_index_to_index(adj_faces[order[0]]);
-    flipped = adj_faces[order[0]] > 0;
+    if (V(F(f, 0),0) == V(F(f, 1),0) && V(F(f, 0),0) == V(F(f, 2),0)) {
+        // The face is perpendicular to X axis.
+        typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+        typedef K::Point_2 Point_2;
+        Point_2 p0(V(F(f,0),1), V(F(f,0),2));
+        Point_2 p1(V(F(f,1),1), V(F(f,1),2));
+        Point_2 p2(V(F(f,2),1), V(F(f,2),2));
+        auto orientation = CGAL::orientation(p0, p1, p2);
+        assert(orientation != CGAL::COLLINEAR);
+        flipped = (orientation == CGAL::NEGATIVE);
+
+    } else {
+        flipped = adj_faces[order[0]] > 0;
+    }
 }