1 // This file is part of libigl, a simple c++ geometry processing library.
2 //
3 // Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public License
6 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
7 // obtain one at http://mozilla.org/MPL/2.0/.
8 #include "order_facets_around_edges.h"
9 #include "order_facets_around_edge.h"
10 #include "../../sort_angles.h"
11 #include <Eigen/Geometry>
12 #include <type_traits>
13 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
14 
15 template<
16     typename DerivedV,
17     typename DerivedF,
18     typename DerivedN,
19     typename DeriveduE,
20     typename uE2EType,
21     typename uE2oEType,
22     typename uE2CType >
23 IGL_INLINE
24 typename std::enable_if<!std::is_same<typename DerivedV::Scalar,
25 typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
order_facets_around_edges(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedN> & N,const Eigen::PlainObjectBase<DeriveduE> & uE,const std::vector<std::vector<uE2EType>> & uE2E,std::vector<std::vector<uE2oEType>> & uE2oE,std::vector<std::vector<uE2CType>> & uE2C)26 igl::copyleft::cgal::order_facets_around_edges(
27         const Eigen::PlainObjectBase<DerivedV>& V,
28         const Eigen::PlainObjectBase<DerivedF>& F,
29         const Eigen::PlainObjectBase<DerivedN>& N,
30         const Eigen::PlainObjectBase<DeriveduE>& uE,
31         const std::vector<std::vector<uE2EType> >& uE2E,
32         std::vector<std::vector<uE2oEType> >& uE2oE,
33         std::vector<std::vector<uE2CType > >& uE2C ) {
34 
35     typedef Eigen::Matrix<typename DerivedN::Scalar, 3, 1> Vector3F;
36     const typename DerivedV::Scalar EPS = 1e-12;
37     const size_t num_faces = F.rows();
38     const size_t num_undirected_edges = uE.rows();
39 
40     auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
41     auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
42 
43     uE2oE.resize(num_undirected_edges);
44     uE2C.resize(num_undirected_edges);
45 
46     for(size_t ui = 0;ui<num_undirected_edges;ui++)
47     {
48         const auto& adj_edges = uE2E[ui];
49         const size_t edge_valance = adj_edges.size();
50         assert(edge_valance > 0);
51 
52         const auto ref_edge = adj_edges[0];
53         const auto ref_face = edge_index_to_face_index(ref_edge);
54         Vector3F ref_normal = N.row(ref_face);
55 
56         const auto ref_corner_o = edge_index_to_corner_index(ref_edge);
57         const auto ref_corner_s = (ref_corner_o+1)%3;
58         const auto ref_corner_d = (ref_corner_o+2)%3;
59 
60         const typename DerivedF::Scalar o = F(ref_face, ref_corner_o);
61         const typename DerivedF::Scalar s = F(ref_face, ref_corner_s);
62         const typename DerivedF::Scalar d = F(ref_face, ref_corner_d);
63 
64         Vector3F edge = V.row(d) - V.row(s);
65         auto edge_len = edge.norm();
66         bool degenerated = edge_len < EPS;
67         if (degenerated) {
68             if (edge_valance <= 2) {
69                 // There is only one way to order 2 or less faces.
70                 edge.setZero();
71             } else {
72                 edge.setZero();
73                 Eigen::Matrix<typename DerivedN::Scalar, Eigen::Dynamic, 3>
74                     normals(edge_valance, 3);
75                 for (size_t fei=0; fei<edge_valance; fei++) {
76                     const auto fe = adj_edges[fei];
77                     const auto f = edge_index_to_face_index(fe);
78                     normals.row(fei) = N.row(f);
79                 }
80                 for (size_t i=0; i<edge_valance; i++) {
81                     size_t j = (i+1) % edge_valance;
82                     Vector3F ni = normals.row(i);
83                     Vector3F nj = normals.row(j);
84                     edge = ni.cross(nj);
85                     edge_len = edge.norm();
86                     if (edge_len >= EPS) {
87                         edge.normalize();
88                         break;
89                     }
90                 }
91 
92                 // Ensure edge direction are consistent with reference face.
93                 Vector3F in_face_vec = V.row(o) - V.row(s);
94                 if (edge.cross(in_face_vec).dot(ref_normal) < 0) {
95                     edge *= -1;
96                 }
97 
98                 if (edge.norm() < EPS) {
99                     std::cerr << "=====================================" << std::endl;
100                     std::cerr << "  ui: " << ui << std::endl;
101                     std::cerr << "edge: " << ref_edge << std::endl;
102                     std::cerr << "face: " << ref_face << std::endl;
103                     std::cerr << "  vs: " << V.row(s) << std::endl;
104                     std::cerr << "  vd: " << V.row(d) << std::endl;
105                     std::cerr << "adj face normals: " << std::endl;
106                     std::cerr << normals << std::endl;
107                     std::cerr << "Very degenerated case detected:" << std::endl;
108                     std::cerr << "Near zero edge surrounded by "
109                         << edge_valance << " neearly colinear faces" <<
110                         std::endl;
111                     std::cerr << "=====================================" << std::endl;
112                 }
113             }
114         } else {
115             edge.normalize();
116         }
117 
118         Eigen::MatrixXd angle_data(edge_valance, 3);
119         std::vector<bool> cons(edge_valance);
120 
121         for (size_t fei=0; fei<edge_valance; fei++) {
122             const auto fe = adj_edges[fei];
123             const auto f = edge_index_to_face_index(fe);
124             const auto c = edge_index_to_corner_index(fe);
125             cons[fei] = (d == F(f, (c+1)%3));
126             assert( cons[fei] ||  (d == F(f,(c+2)%3)));
127             assert(!cons[fei] || (s == F(f,(c+2)%3)));
128             assert(!cons[fei] || (d == F(f,(c+1)%3)));
129             Vector3F n = N.row(f);
130             angle_data(fei, 0) = ref_normal.cross(n).dot(edge);
131             angle_data(fei, 1) = ref_normal.dot(n);
132             if (cons[fei]) {
133                 angle_data(fei, 0) *= -1;
134                 angle_data(fei, 1) *= -1;
135             }
136             angle_data(fei, 0) *= -1; // Sort clockwise.
137             angle_data(fei, 2) = (cons[fei]?1.:-1.)*(f+1);
138         }
139 
140         Eigen::VectorXi order;
141         igl::sort_angles(angle_data, order);
142 
143         auto& ordered_edges = uE2oE[ui];
144         auto& consistency = uE2C[ui];
145 
146         ordered_edges.resize(edge_valance);
147         consistency.resize(edge_valance);
148         for (size_t fei=0; fei<edge_valance; fei++) {
149             ordered_edges[fei] = adj_edges[order[fei]];
150             consistency[fei] = cons[order[fei]];
151         }
152     }
153 }
154 
155 template<
156     typename DerivedV,
157     typename DerivedF,
158     typename DerivedN,
159     typename DeriveduE,
160     typename uE2EType,
161     typename uE2oEType,
162     typename uE2CType >
163 IGL_INLINE
164 typename std::enable_if<std::is_same<typename DerivedV::Scalar,
165 typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
order_facets_around_edges(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedN> & N,const Eigen::PlainObjectBase<DeriveduE> & uE,const std::vector<std::vector<uE2EType>> & uE2E,std::vector<std::vector<uE2oEType>> & uE2oE,std::vector<std::vector<uE2CType>> & uE2C)166 igl::copyleft::cgal::order_facets_around_edges(
167         const Eigen::PlainObjectBase<DerivedV>& V,
168         const Eigen::PlainObjectBase<DerivedF>& F,
169         const Eigen::PlainObjectBase<DerivedN>& N,
170         const Eigen::PlainObjectBase<DeriveduE>& uE,
171         const std::vector<std::vector<uE2EType> >& uE2E,
172         std::vector<std::vector<uE2oEType> >& uE2oE,
173         std::vector<std::vector<uE2CType > >& uE2C ) {
174 
175     typedef Eigen::Matrix<typename DerivedN::Scalar, 3, 1> Vector3F;
176     typedef Eigen::Matrix<typename DerivedV::Scalar, 3, 1> Vector3E;
177     const typename DerivedV::Scalar EPS = 1e-12;
178     const size_t num_faces = F.rows();
179     const size_t num_undirected_edges = uE.rows();
180 
181     auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
182     auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
183 
184     uE2oE.resize(num_undirected_edges);
185     uE2C.resize(num_undirected_edges);
186 
187     for(size_t ui = 0;ui<num_undirected_edges;ui++)
188     {
189         const auto& adj_edges = uE2E[ui];
190         const size_t edge_valance = adj_edges.size();
191         assert(edge_valance > 0);
192 
193         const auto ref_edge = adj_edges[0];
194         const auto ref_face = edge_index_to_face_index(ref_edge);
195         Vector3F ref_normal = N.row(ref_face);
196 
197         const auto ref_corner_o = edge_index_to_corner_index(ref_edge);
198         const auto ref_corner_s = (ref_corner_o+1)%3;
199         const auto ref_corner_d = (ref_corner_o+2)%3;
200 
201         const typename DerivedF::Scalar o = F(ref_face, ref_corner_o);
202         const typename DerivedF::Scalar s = F(ref_face, ref_corner_s);
203         const typename DerivedF::Scalar d = F(ref_face, ref_corner_d);
204 
205         Vector3E exact_edge = V.row(d) - V.row(s);
206         exact_edge.array() /= exact_edge.squaredNorm();
207         Vector3F edge(
208                 CGAL::to_double(exact_edge[0]),
209                 CGAL::to_double(exact_edge[1]),
210                 CGAL::to_double(exact_edge[2]));
211         edge.normalize();
212 
213         Eigen::MatrixXd angle_data(edge_valance, 3);
214         std::vector<bool> cons(edge_valance);
215 
216         for (size_t fei=0; fei<edge_valance; fei++) {
217             const auto fe = adj_edges[fei];
218             const auto f = edge_index_to_face_index(fe);
219             const auto c = edge_index_to_corner_index(fe);
220             cons[fei] = (d == F(f, (c+1)%3));
221             assert( cons[fei] ||  (d == F(f,(c+2)%3)));
222             assert(!cons[fei] || (s == F(f,(c+2)%3)));
223             assert(!cons[fei] || (d == F(f,(c+1)%3)));
224             Vector3F n = N.row(f);
225             angle_data(fei, 0) = ref_normal.cross(n).dot(edge);
226             angle_data(fei, 1) = ref_normal.dot(n);
227             if (cons[fei]) {
228                 angle_data(fei, 0) *= -1;
229                 angle_data(fei, 1) *= -1;
230             }
231             angle_data(fei, 0) *= -1; // Sort clockwise.
232             angle_data(fei, 2) = (cons[fei]?1.:-1.)*(f+1);
233         }
234 
235         Eigen::VectorXi order;
236         igl::sort_angles(angle_data, order);
237 
238         auto& ordered_edges = uE2oE[ui];
239         auto& consistency = uE2C[ui];
240 
241         ordered_edges.resize(edge_valance);
242         consistency.resize(edge_valance);
243         for (size_t fei=0; fei<edge_valance; fei++) {
244             ordered_edges[fei] = adj_edges[order[fei]];
245             consistency[fei] = cons[order[fei]];
246         }
247     }
248 }
249 
250 template<
251     typename DerivedV,
252     typename DerivedF,
253     typename DeriveduE,
254     typename uE2EType,
255     typename uE2oEType,
256     typename uE2CType >
order_facets_around_edges(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DeriveduE> & uE,const std::vector<std::vector<uE2EType>> & uE2E,std::vector<std::vector<uE2oEType>> & uE2oE,std::vector<std::vector<uE2CType>> & uE2C)257 IGL_INLINE void igl::copyleft::cgal::order_facets_around_edges(
258         const Eigen::PlainObjectBase<DerivedV>& V,
259         const Eigen::PlainObjectBase<DerivedF>& F,
260         const Eigen::PlainObjectBase<DeriveduE>& uE,
261         const std::vector<std::vector<uE2EType> >& uE2E,
262         std::vector<std::vector<uE2oEType> >& uE2oE,
263         std::vector<std::vector<uE2CType > >& uE2C ) {
264 
265     //typedef Eigen::Matrix<typename DerivedV::Scalar, 3, 1> Vector3E;
266     const size_t num_faces = F.rows();
267     const size_t num_undirected_edges = uE.rows();
268 
269     auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
270     auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
271 
272     uE2oE.resize(num_undirected_edges);
273     uE2C.resize(num_undirected_edges);
274 
275     for(size_t ui = 0;ui<num_undirected_edges;ui++)
276     {
277         const auto& adj_edges = uE2E[ui];
278         const size_t edge_valance = adj_edges.size();
279         assert(edge_valance > 0);
280 
281         const auto ref_edge = adj_edges[0];
282         const auto ref_face = edge_index_to_face_index(ref_edge);
283 
284         const auto ref_corner_o = edge_index_to_corner_index(ref_edge);
285         const auto ref_corner_s = (ref_corner_o+1)%3;
286         const auto ref_corner_d = (ref_corner_o+2)%3;
287 
288         //const typename DerivedF::Scalar o = F(ref_face, ref_corner_o);
289         const typename DerivedF::Scalar s = F(ref_face, ref_corner_s);
290         const typename DerivedF::Scalar d = F(ref_face, ref_corner_d);
291 
292         std::vector<bool> cons(edge_valance);
293         std::vector<int> adj_faces(edge_valance);
294         for (size_t fei=0; fei<edge_valance; fei++) {
295             const auto fe = adj_edges[fei];
296             const auto f = edge_index_to_face_index(fe);
297             const auto c = edge_index_to_corner_index(fe);
298             cons[fei] = (d == F(f, (c+1)%3));
299             adj_faces[fei] = (f+1) * (cons[fei] ? 1:-1);
300 
301             assert( cons[fei] ||  (d == F(f,(c+2)%3)));
302             assert(!cons[fei] || (s == F(f,(c+2)%3)));
303             assert(!cons[fei] || (d == F(f,(c+1)%3)));
304         }
305 
306         Eigen::VectorXi order;
307         order_facets_around_edge(V, F, s, d, adj_faces, order);
308         assert((size_t)order.size() == edge_valance);
309 
310         auto& ordered_edges = uE2oE[ui];
311         auto& consistency = uE2C[ui];
312 
313         ordered_edges.resize(edge_valance);
314         consistency.resize(edge_valance);
315         for (size_t fei=0; fei<edge_valance; fei++) {
316             ordered_edges[fei] = adj_edges[order[fei]];
317             consistency[fei] = cons[order[fei]];
318         }
319     }
320 }
321 
322 #ifdef IGL_STATIC_LIBRARY
323 // Explicit template instantiation
324 // generated by autoexplicit.sh
325 template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, int, bool>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
326 // generated by autoexplicit.sh
327 template std::enable_if<!(std::is_same<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, CGAL::Lazy_exact_nt<CGAL::Gmpq> >::value), void>::type igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, int, bool>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
328 template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
329 template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
330 template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
331 template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
332 #endif
333