1 // This file is part of libigl, a simple c++ geometry processing library.
2 //
3 // Copyright (C) 2015 Qingnan Zhou <qnzhou@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 "points_inside_component.h"
9 #include "../../LinSpaced.h"
10 #include "order_facets_around_edge.h"
11 #include "assign_scalar.h"
12
13 #include <CGAL/AABB_tree.h>
14 #include <CGAL/AABB_traits.h>
15 #include <CGAL/AABB_triangle_primitive.h>
16 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
17
18 #include <cassert>
19 #include <list>
20 #include <limits>
21 #include <vector>
22
23
24 namespace igl {
25 namespace copyleft
26 {
27 namespace cgal {
28 namespace points_inside_component_helper {
29 typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
30 typedef Kernel::Ray_3 Ray_3;
31 typedef Kernel::Point_3 Point_3;
32 typedef Kernel::Vector_3 Vector_3;
33 typedef Kernel::Triangle_3 Triangle;
34 typedef Kernel::Plane_3 Plane_3;
35 typedef std::vector<Triangle>::iterator Iterator;
36 typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
37 typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
38 typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
39
40 template<typename DerivedF, typename DerivedI>
extract_adj_faces(const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedI> & I,const size_t s,const size_t d,std::vector<int> & adj_faces)41 void extract_adj_faces(
42 const Eigen::PlainObjectBase<DerivedF>& F,
43 const Eigen::PlainObjectBase<DerivedI>& I,
44 const size_t s, const size_t d,
45 std::vector<int>& adj_faces) {
46 const size_t num_faces = I.rows();
47 for (size_t i=0; i<num_faces; i++) {
48 Eigen::Vector3i f = F.row(I(i, 0));
49 if (((size_t)f[0] == s && (size_t)f[1] == d) ||
50 ((size_t)f[1] == s && (size_t)f[2] == d) ||
51 ((size_t)f[2] == s && (size_t)f[0] == d)) {
52 adj_faces.push_back((I(i, 0)+1) * -1);
53 continue;
54 }
55 if (((size_t)f[0] == d && (size_t)f[1] == s) ||
56 ((size_t)f[1] == d && (size_t)f[2] == s) ||
57 ((size_t)f[2] == d && (size_t)f[0] == s)) {
58 adj_faces.push_back(I(i, 0)+1);
59 continue;
60 }
61 }
62 }
63
64 template<typename DerivedF, typename DerivedI>
extract_adj_vertices(const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedI> & I,const size_t v,std::vector<int> & adj_vertices)65 void extract_adj_vertices(
66 const Eigen::PlainObjectBase<DerivedF>& F,
67 const Eigen::PlainObjectBase<DerivedI>& I,
68 const size_t v, std::vector<int>& adj_vertices) {
69 std::set<size_t> unique_adj_vertices;
70 const size_t num_faces = I.rows();
71 for (size_t i=0; i<num_faces; i++) {
72 Eigen::Vector3i f = F.row(I(i, 0));
73 if ((size_t)f[0] == v) {
74 unique_adj_vertices.insert(f[1]);
75 unique_adj_vertices.insert(f[2]);
76 } else if ((size_t)f[1] == v) {
77 unique_adj_vertices.insert(f[0]);
78 unique_adj_vertices.insert(f[2]);
79 } else if ((size_t)f[2] == v) {
80 unique_adj_vertices.insert(f[0]);
81 unique_adj_vertices.insert(f[1]);
82 }
83 }
84 adj_vertices.resize(unique_adj_vertices.size());
85 std::copy(unique_adj_vertices.begin(),
86 unique_adj_vertices.end(),
87 adj_vertices.begin());
88 }
89
90 template<typename DerivedV, typename DerivedF, typename DerivedI>
determine_point_edge_orientation(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedI> & I,const Point_3 & query,size_t s,size_t d)91 bool determine_point_edge_orientation(
92 const Eigen::PlainObjectBase<DerivedV>& V,
93 const Eigen::PlainObjectBase<DerivedF>& F,
94 const Eigen::PlainObjectBase<DerivedI>& I,
95 const Point_3& query, size_t s, size_t d) {
96 // Algorithm:
97 //
98 // Order the adj faces around the edge (s,d) clockwise using
99 // query point as pivot. (i.e. The first face of the ordering
100 // is directly after the pivot point, and the last face is
101 // directly before the pivot.)
102 //
103 // The point is outside if the first and last faces of the
104 // ordering forms a convex angle. This check can be done
105 // without any construction by looking at the orientation of the
106 // faces. The angle is convex iff the first face contains (s,d)
107 // as an edge and the last face contains (d,s) as an edge.
108 //
109 // The point is inside if the first and last faces of the
110 // ordering forms a concave angle. That is the first face
111 // contains (d,s) as an edge and the last face contains (s,d) as
112 // an edge.
113 //
114 // In the special case of duplicated faces. I.e. multiple faces
115 // sharing the same 3 corners, but not necessarily the same
116 // orientation. The ordering will always rank faces containing
117 // edge (s,d) before faces containing edge (d,s).
118 //
119 // Therefore, if there are any duplicates of the first faces,
120 // the ordering will always choose the one with edge (s,d) if
121 // possible. The same for the last face.
122 //
123 // In the very degenerated case where the first and last face
124 // are duplicates, but with different orientations, it is
125 // equally valid to think the angle formed by them is either 0
126 // or 360 degrees. By default, 0 degree is used, and thus the
127 // query point is outside.
128
129 std::vector<int> adj_faces;
130 extract_adj_faces(F, I, s, d, adj_faces);
131 const size_t num_adj_faces = adj_faces.size();
132 assert(num_adj_faces > 0);
133
134 DerivedV pivot_point(1, 3);
135 igl::copyleft::cgal::assign_scalar(query.x(), pivot_point(0, 0));
136 igl::copyleft::cgal::assign_scalar(query.y(), pivot_point(0, 1));
137 igl::copyleft::cgal::assign_scalar(query.z(), pivot_point(0, 2));
138 Eigen::VectorXi order;
139 order_facets_around_edge(V, F, s, d,
140 adj_faces, pivot_point, order);
141 assert((size_t)order.size() == num_adj_faces);
142 if (adj_faces[order[0]] > 0 &&
143 adj_faces[order[num_adj_faces-1] < 0]) {
144 return true;
145 } else if (adj_faces[order[0]] < 0 &&
146 adj_faces[order[num_adj_faces-1] > 0]) {
147 return false;
148 } else {
149 throw "The input mesh does not represent a valid volume";
150 }
151 throw "The input mesh does not represent a valid volume";
152 return false;
153 }
154
155 template<typename DerivedV, typename DerivedF, typename DerivedI>
determine_point_vertex_orientation(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedI> & I,const Point_3 & query,size_t s)156 bool determine_point_vertex_orientation(
157 const Eigen::PlainObjectBase<DerivedV>& V,
158 const Eigen::PlainObjectBase<DerivedF>& F,
159 const Eigen::PlainObjectBase<DerivedI>& I,
160 const Point_3& query, size_t s) {
161 std::vector<int> adj_vertices;
162 extract_adj_vertices(F, I, s, adj_vertices);
163 const size_t num_adj_vertices = adj_vertices.size();
164
165 std::vector<Point_3> adj_points;
166 for (size_t i=0; i<num_adj_vertices; i++) {
167 const size_t vi = adj_vertices[i];
168 adj_points.emplace_back(V(vi,0), V(vi,1), V(vi,2));
169 }
170
171 // A plane is on the exterior if all adj_points lies on or to
172 // one side of the plane.
173 auto is_on_exterior = [&](const Plane_3& separator) -> bool{
174 size_t positive=0;
175 size_t negative=0;
176 size_t coplanar=0;
177 for (const auto& point : adj_points) {
178 switch(separator.oriented_side(point)) {
179 case CGAL::ON_POSITIVE_SIDE:
180 positive++;
181 break;
182 case CGAL::ON_NEGATIVE_SIDE:
183 negative++;
184 break;
185 case CGAL::ON_ORIENTED_BOUNDARY:
186 coplanar++;
187 break;
188 default:
189 throw "Unknown plane-point orientation";
190 }
191 }
192 auto query_orientation = separator.oriented_side(query);
193 bool r =
194 (positive == 0 && query_orientation == CGAL::POSITIVE)
195 ||
196 (negative == 0 && query_orientation == CGAL::NEGATIVE);
197 return r;
198 };
199
200 size_t d = std::numeric_limits<size_t>::max();
201 Point_3 p(V(s,0), V(s,1), V(s,2));
202 for (size_t i=0; i<num_adj_vertices; i++) {
203 const size_t vi = adj_vertices[i];
204 for (size_t j=i+1; j<num_adj_vertices; j++) {
205 Plane_3 separator(p, adj_points[i], adj_points[j]);
206 if (separator.is_degenerate()) {
207 throw "Input mesh contains degenerated faces";
208 }
209 if (is_on_exterior(separator)) {
210 d = vi;
211 assert(!CGAL::collinear(p, adj_points[i], query));
212 break;
213 }
214 }
215 if (d < (size_t)V.rows()) break;
216 }
217 if (d > (size_t)V.rows()) {
218 // All adj faces are coplanar, use the first edge.
219 d = adj_vertices[0];
220 }
221 return determine_point_edge_orientation(V, F, I, query, s, d);
222 }
223
224 template<typename DerivedV, typename DerivedF, typename DerivedI>
determine_point_face_orientation(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedI> & I,const Point_3 & query,size_t fid)225 bool determine_point_face_orientation(
226 const Eigen::PlainObjectBase<DerivedV>& V,
227 const Eigen::PlainObjectBase<DerivedF>& F,
228 const Eigen::PlainObjectBase<DerivedI>& I,
229 const Point_3& query, size_t fid) {
230 // Algorithm: A point is on the inside of a face if the
231 // tetrahedron formed by them is negatively oriented.
232
233 Eigen::Vector3i f = F.row(I(fid, 0));
234 const Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
235 const Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
236 const Point_3 v2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
237 auto result = CGAL::orientation(v0, v1, v2, query);
238 if (result == CGAL::COPLANAR) {
239 throw "Cannot determine inside/outside because query point lies exactly on the input surface.";
240 }
241 return result == CGAL::NEGATIVE;
242 }
243 }
244 }
245 }
246 }
247
248 template<typename DerivedV, typename DerivedF, typename DerivedI,
249 typename DerivedP, typename DerivedB>
points_inside_component(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedI> & I,const Eigen::PlainObjectBase<DerivedP> & P,Eigen::PlainObjectBase<DerivedB> & inside)250 IGL_INLINE void igl::copyleft::cgal::points_inside_component(
251 const Eigen::PlainObjectBase<DerivedV>& V,
252 const Eigen::PlainObjectBase<DerivedF>& F,
253 const Eigen::PlainObjectBase<DerivedI>& I,
254 const Eigen::PlainObjectBase<DerivedP>& P,
255 Eigen::PlainObjectBase<DerivedB>& inside) {
256 using namespace igl::copyleft::cgal::points_inside_component_helper;
257 if (F.rows() <= 0 || I.rows() <= 0) {
258 throw "Inside check cannot be done on empty facet component.";
259 }
260
261 const size_t num_faces = I.rows();
262 std::vector<Triangle> triangles;
263 for (size_t i=0; i<num_faces; i++) {
264 const Eigen::Vector3i f = F.row(I(i, 0));
265 triangles.emplace_back(
266 Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
267 Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
268 Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
269 if (triangles.back().is_degenerate()) {
270 throw "Input facet components contains degenerated triangles";
271 }
272 }
273 Tree tree(triangles.begin(), triangles.end());
274 tree.accelerate_distance_queries();
275
276 enum ElementType { VERTEX, EDGE, FACE };
277 auto determine_element_type = [&](
278 size_t fid, const Point_3& p, size_t& element_index) -> ElementType{
279 const Eigen::Vector3i f = F.row(I(fid, 0));
280 const Point_3 p0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
281 const Point_3 p1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
282 const Point_3 p2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
283
284 if (p == p0) { element_index = 0; return VERTEX; }
285 if (p == p1) { element_index = 1; return VERTEX; }
286 if (p == p2) { element_index = 2; return VERTEX; }
287 if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
288 if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
289 if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
290
291 element_index = 0;
292 return FACE;
293 };
294
295 const size_t num_queries = P.rows();
296 inside.resize(num_queries, 1);
297 for (size_t i=0; i<num_queries; i++) {
298 const Point_3 query(P(i,0), P(i,1), P(i,2));
299 auto projection = tree.closest_point_and_primitive(query);
300 auto closest_point = projection.first;
301 size_t fid = projection.second - triangles.begin();
302
303 size_t element_index;
304 switch (determine_element_type(fid, closest_point, element_index)) {
305 case VERTEX:
306 {
307 const size_t s = F(I(fid, 0), element_index);
308 inside(i,0) = determine_point_vertex_orientation(
309 V, F, I, query, s);
310 }
311 break;
312 case EDGE:
313 {
314 const size_t s = F(I(fid, 0), (element_index+1)%3);
315 const size_t d = F(I(fid, 0), (element_index+2)%3);
316 inside(i,0) = determine_point_edge_orientation(
317 V, F, I, query, s, d);
318 }
319 break;
320 case FACE:
321 inside(i,0) = determine_point_face_orientation(V, F, I, query, fid);
322 break;
323 default:
324 throw "Unknown closest element type!";
325 }
326 }
327 }
328
329 template<typename DerivedV, typename DerivedF, typename DerivedP,
330 typename DerivedB>
points_inside_component(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedP> & P,Eigen::PlainObjectBase<DerivedB> & inside)331 IGL_INLINE void igl::copyleft::cgal::points_inside_component(
332 const Eigen::PlainObjectBase<DerivedV>& V,
333 const Eigen::PlainObjectBase<DerivedF>& F,
334 const Eigen::PlainObjectBase<DerivedP>& P,
335 Eigen::PlainObjectBase<DerivedB>& inside) {
336 Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi>(F.rows(), 0, F.rows()-1);
337 igl::copyleft::cgal::points_inside_component(V, F, I, P, inside);
338 }
339
340 #ifdef IGL_STATIC_LIBRARY
341 // Explicit template instantiation
342 // generated by autoexplicit.sh
343 template void igl::copyleft::cgal::points_inside_component<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Array<bool, -1, 1, 0, -1, 1> >(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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Array<bool, -1, 1, 0, -1, 1> >&);
344 template void igl::copyleft::cgal::points_inside_component< Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix< int, -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> > ( 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&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix< int, -1, -1, 0, -1, -1> >&);
345 template void igl::copyleft::cgal::points_inside_component< 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> > ( 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> >&);
346 template void igl::copyleft::cgal::points_inside_component<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
347 template void igl::copyleft::cgal::points_inside_component<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(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<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
348 template void igl::copyleft::cgal::points_inside_component<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> >(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> >&);
349 template void igl::copyleft::cgal::points_inside_component<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
350 #endif
351