Added some commented out code.

Former-commit-id: 0ddbf149fe4c1b85d01b1abf3533626a65b8c300

include/igl/cgal/propagate_winding_numbers.cpp

@@ -3,9 +3,12 @@
 #include "../extract_non_manifold_edge_curves.h"
 #include "../facet_components.h"
 #include "../unique_edge_map.h"
+#include "../writeOBJ.h"
 #include "order_facets_around_edge.h"
 #include "outer_facet.h"
 #include "closest_facet.h"
+#include "assign_scalar.h"
+#include "extract_cells.h"
 
 #include <stdexcept>
 #include <limits>
@@ -38,6 +41,9 @@ namespace propagate_winding_numbers_helper {
                     int next_winding_number = per_patch_winding_number(
                             order[next], k*2+(orientation[next]? 1:0));
                     if (curr_winding_number != next_winding_number) {
+                        std::cout << "edge: " << i << std::endl;
+                        std::cout << curr_winding_number << " != " <<
+                            next_winding_number << std::endl;
                         return false;
                     }
                 }
@@ -185,6 +191,14 @@ IGL_INLINE bool igl::cgal::propagate_winding_numbers_single_component_patch_wise
             assert(curr_i < num_adj_patches);
             const bool curr_ori = orientation[curr_i];
 
+            //for (size_t i=0; i<num_adj_patches; i++) {
+            //    const size_t next_i = (curr_i + 1) % num_adj_patches;
+            //    const size_t num_patch_idx = order[next_i];
+            //    const bool next_ori = orientation[next_i];
+            //    const size_t next_patch_label = labels[next_patch_idx];
+            //    // TODO
+            //}
+
             const size_t next_i = curr_ori ? (curr_i+1) % num_adj_patches
                 : (curr_i+num_adj_patches-1)%num_adj_patches;
             const size_t prev_i = curr_ori ?
@@ -250,10 +264,298 @@ IGL_INLINE bool igl::cgal::propagate_winding_numbers_single_component_patch_wise
             }
         }
     }
+    assert((patch_W.array() != INVALID).all());
+
+    using namespace propagate_winding_numbers_helper;
+    return winding_number_assignment_is_consistent(orders, orientations, patch_W);
+}
+
+#if 0
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DeriveduE,
+    typename uE2EType,
+    typename DerivedEMAP,
+    typename DerivedC,
+    typename DerivedP,
+    typename DerivedW >
+IGL_INLINE bool igl::cgal::propagate_winding_numbers_single_component_patch_wise(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        const Eigen::PlainObjectBase<DerivedC>& labels,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedW>& patch_W) {
+    const size_t num_faces = F.rows();
+    const size_t num_patches = P.maxCoeff() + 1;
+    assert(labels.size() == num_patches);
+    // Utility functions.
+    auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
+    auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
+    auto face_and_corner_index_to_edge_index = [&](size_t fi, size_t ci) {
+        return ci*num_faces + fi;
+    };
+    auto get_edge_end_points = [&](size_t ei, size_t& s, size_t& d) {
+        const size_t fi = edge_index_to_face_index(ei);
+        const size_t ci = edge_index_to_corner_index(ei);
+        s = F(fi, (ci+1)%3);
+        d = F(fi, (ci+2)%3);
+    };
+    auto is_positively_orientated =
+        [&](size_t fi, size_t s, size_t d) -> bool{
+            const Eigen::Vector3i f = F.row(fi);
+            if (f[0] == d && f[1] == s) return true;
+            if (f[1] == d && f[2] == s) return true;
+            if (f[2] == d && f[0] == s) return true;
+            if (f[0] == s && f[1] == d) return false;
+            if (f[1] == s && f[2] == d) return false;
+            if (f[2] == s && f[0] == d) return false;
+            throw std::runtime_error("Edge does not belong to face.");
+            return false;
+        };
+
+    auto compute_signed_index = [&](size_t fi, size_t s, size_t d) -> int{
+        return int(fi+1) * (is_positively_orientated(fi, s, d) ? 1:-1);
+    };
+    auto compute_unsigned_index = [&](int signed_idx) -> size_t {
+        return abs(signed_idx) - 1;
+    };
+
+    // Order patches around each non-manifold edge.
+    const size_t num_edges = uE.rows();
+    std::vector<Eigen::VectorXi> orders(num_edges);
+    std::vector<std::vector<bool> > orientations(num_edges);
+    std::vector<std::vector<size_t> > patch_edge_adjacency(num_patches);
+    for (size_t uei=0; uei<num_edges; uei++) {
+        const auto& edges = uE2E[uei];
+        if (edges.size() <= 2) continue;
+        size_t s = uE(uei, 0);
+        size_t d = uE(uei, 1);
+        std::vector<int> adj_faces;
+        for (size_t ei : uE2E[uei]) {
+            const size_t fi = edge_index_to_face_index(ei);
+            const size_t signed_fi =
+                compute_signed_index(fi, s, d);
+            adj_faces.push_back(signed_fi);
+            patch_edge_adjacency[P[fi]].push_back(ei);
+        }
+
+        auto& order = orders[uei];
+        igl::cgal::order_facets_around_edge(
+                V, F, s, d, adj_faces, order);
+        assert(order.minCoeff() == 0);
+        assert(order.maxCoeff() == adj_faces.size() - 1);
+
+        auto& orientation = orientations[uei];
+        orientation.resize(order.size());
+        std::transform(order.data(), order.data()+order.size(),
+                orientation.begin(), 
+                [&](size_t index) { return adj_faces[index] > 0; });
+        std::transform(order.data(), order.data()+order.size(),
+                order.data(), [&](size_t index) {
+                return compute_unsigned_index(adj_faces[index]);
+                } );
+    }
+
+    // Propagate winding number from infinity.
+    // Assuming infinity has winding number 0.
+    const size_t num_labels = labels.maxCoeff() + 1;
+    const int INVALID = std::numeric_limits<int>::max();
+    patch_W.resize(num_patches, 2*num_labels);
+    patch_W.setConstant(INVALID);
+
+    size_t outer_facet_idx;
+    bool outer_facet_is_flipped;
+    Eigen::VectorXi face_indices(num_faces);
+    face_indices.setLinSpaced(num_faces, 0, num_faces-1);
+    igl::cgal::outer_facet(V, F, face_indices,
+            outer_facet_idx, outer_facet_is_flipped);
+    size_t outer_patch_idx = P[outer_facet_idx];
+    size_t outer_patch_label = labels[outer_patch_idx];
+    patch_W.row(outer_patch_idx).setZero();
+    if (outer_facet_is_flipped) {
+        patch_W(outer_patch_idx, outer_patch_label*2) = -1;
+    } else {
+        patch_W(outer_patch_idx, outer_patch_label*2+1) = 1;
+    }
+
+    auto winding_num_assigned = [&](size_t patch_idx) -> bool{
+        return (patch_W.row(patch_idx).array() != INVALID).all();
+    };
+
+    std::queue<size_t> Q;
+    Q.push(outer_patch_idx);
+    while (!Q.empty()) {
+        const size_t curr_patch_idx = Q.front();
+        const size_t curr_patch_label = labels[curr_patch_idx];
+        Q.pop();
+
+        const auto& adj_edges = patch_edge_adjacency[curr_patch_idx];
+        for (size_t ei: adj_edges) {
+            const size_t uei = EMAP[ei];
+            const size_t fid = edge_index_to_face_index(ei);
+            const auto& order = orders[uei];
+            const auto& orientation = orientations[uei];
+            const size_t num_adj_patches = order.size();
+            assert(num_adj_patches == orientation.size());
+
+            size_t curr_i = std::numeric_limits<size_t>::max();
+            for (size_t i=0; i<num_adj_patches; i++) {
+                std::cout << (orientation[i]?"+":"-") << order[i] << " ";
+            }
+            std::cout << std::endl;
+            for (size_t i=0; i<num_adj_patches; i++) {
+                if (order[i] == fid) {
+                    curr_i = i;
+                    break;
+                }
+            }
+            assert(curr_i < num_adj_patches);
+            const bool curr_ori = orientation[curr_i];
+
+            const size_t next_i = curr_ori ? (curr_i+1) % num_adj_patches
+                : (curr_i+num_adj_patches-1)%num_adj_patches;
+            const size_t prev_i = curr_ori ?
+                (curr_i+num_adj_patches-1)%num_adj_patches
+                : (curr_i+1) % num_adj_patches;
+            const size_t next_patch_idx = P[order[next_i]];
+            const size_t prev_patch_idx = P[order[prev_i]];
+
+            if (!winding_num_assigned(next_patch_idx)) {
+                const bool next_ori = orientation[next_i];
+                const bool next_cons = next_ori != curr_ori;
+                const size_t next_patch_label = labels[next_patch_idx];
+                for (size_t i=0; i<num_labels; i++) {
+                    const int shared_winding_number =
+                        patch_W(curr_patch_idx, i*2);
+
+                    if (i == next_patch_label) {
+                        // Truth table:
+                        // curr_ori  next_ori  wind_# inc
+                        // True      True      -1
+                        // True      False      1
+                        // False     True       1
+                        // False     False     -1
+
+                        patch_W(next_patch_idx, i*2+(next_cons ?0:1)) =
+                            shared_winding_number;
+                        patch_W(next_patch_idx, i*2+(next_cons ?1:0)) = 
+                            shared_winding_number + (next_cons ? 1:-1);
+                    } else {
+                        patch_W(next_patch_idx, i*2  ) = shared_winding_number;
+                        patch_W(next_patch_idx, i*2+1) = shared_winding_number;
+                    }
+                }
+                Q.push(next_patch_idx);
+            } else {
+                const bool next_ori = orientation[next_i];
+                const bool next_cons = next_ori != curr_ori;
+                const size_t next_patch_label = labels[next_patch_idx];
+                for (size_t i=0; i<num_labels; i++) {
+                    const int shared_winding_number =
+                        patch_W(curr_patch_idx, i*2);
+
+                    if (i == next_patch_label) {
+                        // Truth table:
+                        // curr_ori  next_ori  wind_# inc
+                        // True      True      -1
+                        // True      False      1
+                        // False     True       1
+                        // False     False     -1
+
+                        assert(patch_W(next_patch_idx, i*2+(next_cons ?0:1)) ==
+                            shared_winding_number);
+                        assert(patch_W(next_patch_idx, i*2+(next_cons ?1:0)) == 
+                            shared_winding_number + (next_cons ? 1:-1));
+                    } else {
+                        assert(patch_W(next_patch_idx, i*2) ==
+                                shared_winding_number);
+                        assert(patch_W(next_patch_idx, i*2+1) ==
+                                shared_winding_number);
+                    }
+                }
+            }
+            if (!winding_num_assigned(prev_patch_idx)) {
+                const bool prev_ori = orientation[prev_i];
+                const bool prev_cons = prev_ori != curr_ori;
+                const size_t prev_patch_label = labels[prev_patch_idx];
+
+                for (size_t i=0; i<num_labels; i++) {
+                    const int shared_winding_number =
+                        patch_W(curr_patch_idx, i*2+1);
+
+                    if (i == prev_patch_label) {
+                        // Truth table:
+                        // curr_ori  next_ori  wind_# inc
+                        // True      True       1
+                        // True      False     -1
+                        // False     True      -1
+                        // False     False      1
+
+                        patch_W(prev_patch_idx, i*2+(prev_cons ?1:0)) =
+                            shared_winding_number;
+                        patch_W(prev_patch_idx, i*2+(prev_cons ?0:1)) =
+                            shared_winding_number + (prev_cons ? -1:1);
+                    } else {
+                        patch_W(prev_patch_idx, i*2  ) = shared_winding_number; 
+                        patch_W(prev_patch_idx, i*2+1) = shared_winding_number; 
+                    }
+                }
+                Q.push(prev_patch_idx);
+            } else {
+                const bool prev_ori = orientation[prev_i];
+                const bool prev_cons = prev_ori != curr_ori;
+                const size_t prev_patch_label = labels[prev_patch_idx];
+
+                for (size_t i=0; i<num_labels; i++) {
+                    const int shared_winding_number =
+                        patch_W(curr_patch_idx, i*2+1);
+
+                    if (i == prev_patch_label) {
+                        // Truth table:
+                        // curr_ori  next_ori  wind_# inc
+                        // True      True       1
+                        // True      False     -1
+                        // False     True      -1
+                        // False     False      1
+
+                        assert(patch_W(prev_patch_idx, i*2+(prev_cons ?1:0)) ==
+                            shared_winding_number);
+                        assert(patch_W(prev_patch_idx, i*2+(prev_cons ?0:1)) ==
+                            shared_winding_number + (prev_cons ? -1:1));
+                    } else {
+                        assert(patch_W(prev_patch_idx, i*2  ) ==
+                                shared_winding_number); 
+                        assert(patch_W(prev_patch_idx, i*2+1) ==
+                                shared_winding_number); 
+                    }
+                }
+            }
+        }
+    }
+    assert((patch_W.array() != INVALID).all());
 
     using namespace propagate_winding_numbers_helper;
+    if (!winding_number_assignment_is_consistent(orders, orientations, patch_W)) {
+        Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
+            DerivedV::IsRowMajor> Vd(V.rows(), V.cols());
+        for (size_t j=0; j<V.rows(); j++) {
+            for (size_t k=0; k<V.cols(); k++) {
+                igl::cgal::assign_scalar(V(j,k), Vd(j,k));
+            }
+        }
+        //for (size_t j=0; j<num_faces; j++) {
+        //    std::cout << patch_W(P[j], 1) << std::endl;
+        //}
+        igl::writeOBJ("debug_wn.obj", Vd, F);
+        assert(false);
+    }
     return winding_number_assignment_is_consistent(orders, orientations, patch_W);
 }
+#endif
 
 template<
 typename DerivedV,
@@ -305,9 +607,11 @@ IGL_INLINE bool igl::cgal::propagate_winding_numbers_single_component(
     assert(winding_numbers.cols() == 2 * num_labels);
 
     W.resize(num_faces, 2*num_labels);
+    W.setConstant(INVALID);
     for (size_t i=0; i<num_faces; i++) {
         W.row(i) = winding_numbers.row(P[i]);
     }
+    assert((W.array() != INVALID).any());
 
     return is_consistent;
 }
@@ -404,10 +708,20 @@ IGL_INLINE void igl::cgal::propagate_winding_numbers(
         }
 
         if (!is_consistent) {
+            Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
+                DerivedV::IsRowMajor> Vd(V.rows(), V.cols());
+            for (size_t j=0; j<V.rows(); j++) {
+                for (size_t k=0; k<V.cols(); k++) {
+                    igl::cgal::assign_scalar(V(j,k), Vd(j,k));
+                }
+            }
+            igl::writeOBJ("debug_wn.obj", Vd, comp_faces[i]);
+            std::cout << Ws[i].col(1) << std::endl;
             std::stringstream err_msg;
             err_msg << "Component " << i
                 << " has inconsistant winding number assignment." << std::endl;
-            throw std::runtime_error(err_msg.str());
+            //throw std::runtime_error(err_msg.str());
+            std::cout << err_msg.str() << std::endl;
         }
     }
 
@@ -419,31 +733,33 @@ IGL_INLINE void igl::cgal::propagate_winding_numbers(
     std::vector<Eigen::MatrixXi> nested_Ws = Ws;
     Eigen::MatrixXi ambient_correction(num_components, 2*num_labels);
     ambient_correction.setZero();
-    for (size_t i=0; i<num_components; i++) {
-        DerivedV samples(num_components-1, 3);
-        Eigen::VectorXi is_inside;
-        auto index_without_i = [&](size_t index) {
-            return index < i ? index:index-1;
-        };
-        for (size_t j=0; j<num_components; j++) {
-            if (i == j) continue;
-            samples.row(index_without_i(j)) = sample_component(j);
-        }
-        Eigen::VectorXi fids;
-        Eigen::Matrix<bool, Eigen::Dynamic, 1> orientation;
-        igl::cgal::closest_facet(V, comp_faces[i], samples,
-                fids, orientation);
-
-        const auto& comp = components[i];
-        for (size_t j=0; j<num_components; j++) {
-            if (i == j) continue;
-            const size_t index = index_without_i(j);
-            const size_t fid = fids(index, 0);
-            const bool ori = orientation(index, 0);
-            for (size_t k=0; k<num_labels; k++) {
-                const int correction = W(comp[fid], k*2+(ori?0:1));
-                ambient_correction(j, k*2  ) += correction;
-                ambient_correction(j, k*2+1) += correction;
+    if (num_components > 1) {
+        for (size_t i=0; i<num_components; i++) {
+            DerivedV samples(num_components-1, 3);
+            Eigen::VectorXi is_inside;
+            auto index_without_i = [&](size_t index) {
+                return index < i ? index:index-1;
+            };
+            for (size_t j=0; j<num_components; j++) {
+                if (i == j) continue;
+                samples.row(index_without_i(j)) = sample_component(j);
+            }
+            Eigen::VectorXi fids;
+            Eigen::Matrix<bool, Eigen::Dynamic, 1> orientation;
+            igl::cgal::closest_facet(V, comp_faces[i], samples,
+                    fids, orientation);
+
+            const auto& comp = components[i];
+            for (size_t j=0; j<num_components; j++) {
+                if (i == j) continue;
+                const size_t index = index_without_i(j);
+                const size_t fid = fids(index, 0);
+                const bool ori = orientation(index, 0);
+                for (size_t k=0; k<num_labels; k++) {
+                    const int correction = W(comp[fid], k*2+(ori?0:1));
+                    ambient_correction(j, k*2  ) += correction;
+                    ambient_correction(j, k*2+1) += correction;
+                }
             }
         }
     }

include/igl/cgal/propagate_winding_numbers.h

@@ -64,6 +64,27 @@ namespace igl {
                 const std::vector<std::vector<size_t> >& intersection_curves,
                 Eigen::PlainObjectBase<DerivedW>& patch_W);
 
+#if 0
+        template<
+            typename DerivedV,
+            typename DerivedF,
+            typename DeriveduE,
+            typename uE2EType,
+            typename DerivedEMAP,
+            typename DerivedC,
+            typename DerivedP,
+            typename DerivedW >
+        IGL_INLINE bool propagate_winding_numbers_single_component_patch_wise(
+                const Eigen::PlainObjectBase<DerivedV>& V,
+                const Eigen::PlainObjectBase<DerivedF>& F,
+                const Eigen::PlainObjectBase<DeriveduE>& uE,
+                const std::vector<std::vector<uE2EType> >& uE2E,
+                const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+                const Eigen::PlainObjectBase<DerivedC>& labels,
+                const Eigen::PlainObjectBase<DerivedP>& P,
+                Eigen::PlainObjectBase<DerivedW>& patch_W);
+#endif
+
         // Compute winding number on each side of the face.  The input mesh must
         // be a single edge-connected component.
         //