cgv
convex_polyhedron.h
1 #pragma once
2 
3 #include <map>
4 #include <iterator>
5 #include <cgv/math/fvec.h>
6 #include <cgv/media/axis_aligned_box.h>
7 #include <cgv/media/color.h>
8 
9 #include "../lib_begin.h"
10 
11 namespace cgv {
12  namespace media {
13  namespace mesh {
14 
15 class CGV_API convex_polyhedron_base
16 {
17 public:
19  enum VertexPlaneLocation {
20  VPL_INSIDE = -1,
21  VPL_TOUCH = 0,
22  VPL_OUTSIDE = 1
23  };
25  enum FacePlaneLocation {
26  // valid cases
27  FPL_OUTSIDE = 0, // face is outside (is refined by TOUCH-versions by adding or or-ring together)
28  FPL_INSIDE = 3, // face is inside (is refined by TOUCH-versions by adding or or-ring together)
29 
30  FPL_TOUCH_NOTHING = 0, // face is not touching and is not split but completely outside or inside (v#+|-=0, v#0=0)[e#0=0, e#s=0, #t=0]
31  FPL_TOUCH_VERTEX = 1, // one vertex of face touches plane (v#+|-=0, v#0=1)[e#0=0, e#s=0, #t=0]
32  FPL_TOUCH_EDGE = 2, // one edge of face touches plane (v#+|-=0, v#0=2, e#0=1)[e#s=0, #t=0]
33 
34  FPL_TOUCH_FACE = 6, // whole face touches plane (v#+|-=0, v#0=n)[e#0=0, e#s=0, #t=0]
35 
36  FPL_SPLIT_TWO_VERTICES, // face is split through two not adjacent vertices (#t=2, v#0=2, e#s=0)[e#0=0, v#+>0, v#->0]
37  FPL_SPLIT_VERTEX_EDGE, // face is split through one vertex and one edge (#t=2, v#0=1, e#s=1)[e#0=0, v#+>0, v#->0]
38  FPL_SPLIT_TWO_EDGES, // face is split through two edges (#t=2, v#0=0, e#s=2)[e#0=0, v#+>0, v#->0]
39 
40  // causes of invalidity
41  FPL_TWO_NON_ADJACENT_VERTICES_TOUCH,// two non adjacent vertices touch plane (v#+|-=0, v#0=2, e#0=0)
42  FPL_INCOMLETE_TOUCH_FACE, // more than two and less than all vertices touch plane (v#0>2, v#0<n)
43  FPL_TOUCH_AND_SPLIT, // in addition to a split we have further touch events (#t=2, v#0+e#s>2)
44  FPL_MULTI_SPLIT, // more than one pair of inside->outside, outside->inside transitions found (#t>2)
45  FPL_ODD_TRANSITION_NUMBER, // the computed number of transitions is odd (#t%2=1)
46 
47  FPL_UNCONSIDERED_CASE,
48 
49  FPL_FIRST_INVALID_CASE = FPL_TWO_NON_ADJACENT_VERTICES_TOUCH
50  };
51 
53  typedef std::vector<int> face_type;
54 
56  typedef std::vector<face_type> shell_type;
57 
58 protected:
60  std::vector<shell_type> shells;
61 public:
63  static int to_index(VertexPlaneLocation vertex_location) { return (int)vertex_location + 1; }
65  static bool is_face_location_on_one_side(FacePlaneLocation face_location) { return face_location < FPL_TOUCH_FACE; }
67  static bool is_face_location_split(FacePlaneLocation face_location) { return face_location > FPL_TOUCH_FACE && face_location < FPL_FIRST_INVALID_CASE; }
69  static bool is_valid_face_location(FacePlaneLocation face_location) { return face_location < FPL_FIRST_INVALID_CASE; }
71  static int box_get_face_vertex_index(int fi, int ci);
73  static void change_orientation(face_type& face);
75  static VertexPlaneLocation vertex_location_side(const std::vector<VertexPlaneLocation>& vertex_locations);
77  static VertexPlaneLocation shell_location_side(const std::vector<FacePlaneLocation>& face_locations);
78 
80  static std::pair<int, int> extract_touching_halfedge(const face_type& face, const std::vector<VertexPlaneLocation>& vertex_locations);
82  FacePlaneLocation get_face_plane_location(const face_type& face, const std::vector<VertexPlaneLocation>& vertex_locations) const;
85  size_t get_nr_shells() const { return shells.size(); }
88  const shell_type& shell(int si) const { return shells[si]; }
90  shell_type& shell(int si) { return shells[si]; }
92  virtual int add_shell() { shells.push_back(shell_type()); return (int)shells.size() - 1; }
94  virtual void del_shell(int si) { shells.erase(shells.begin() + si); }
96  virtual void copy_shell(int si_source, int si_target) { shell(si_target) = shell(si_source); }
98  virtual void swap_shells(int si0, int si1) { std::swap(shell(si0), shell(si1)); }
100  size_t get_nr_faces(int si = 0) const { return shell(si).size(); }
102  virtual int add_face(int si = 0) { shell(si).push_back(face_type()); return (int)shell(si).size() - 1; }
104  virtual void del_face(int fi, int si = 0) { shell(si).erase(shell(si).begin() + fi); }
106  const face_type& face(int fi, int si = 0) const { return shell(si)[fi]; }
108  face_type& face(int fi, int si = 0) { return shell(si)[fi]; }
110 
112  virtual size_t get_nr_vertices() const = 0;
115  virtual void del_vertex(int vi) = 0;
117  void remove_redundant_vertices();
119 };
120 
122 
123 template <typename T, cgv::type::int32_type TCDim = 2>
124 class convex_polyhedron : public convex_polyhedron_base
125 {
126 public:
128  typedef T coord_type;
130  const cgv::type::int32_type texcoord_dim = TCDim;
132  typedef typename cgv::math::fvec<T, 3> point_type;
134  typedef typename cgv::math::fvec<T, TCDim> texcoord_type;
136  typedef typename cgv::math::fvec<T, 4> plane_type;
138  typedef typename cgv::media::axis_aligned_box<T, 3> box_type;
140  struct vertex_type
141  {
142  point_type position;
143  texcoord_type texcoord;
144  };
145 protected:
147  std::vector<vertex_type> vertices;
149  std::vector<std::vector<plane_type> > face_planes;
150 public:
153  size_t get_nr_vertices() const { return vertices.size(); }
156  virtual int add_vertex() { vertices.push_back(vertex_type()); return (int)get_nr_vertices() - 1; }
158  void del_vertex(int vi) { vertices.erase(vertices.begin() + vi); }
160  const vertex_type& vertex(int vi) const { return vertices[vi]; }
162  vertex_type& vertex(int vi) { return vertices[vi]; }
164  const point_type& position(int vi) const { return vertices[vi].position; }
166  point_type& position(int vi) { return vertices[vi].position; }
168  const texcoord_type& texcoord(int vi) const { return vertices[vi].texcoord; }
170  texcoord_type& texcoord(int vi) { return vertices[vi].texcoord; }
172 
173 
176  static bool orientation_match(const plane_type& p, const plane_type& q) {
178  return dot(reinterpret_cast<const point_type&>(p), reinterpret_cast<const point_type&>(q)) > 0;
179  }
181  const plane_type& face_plane(int fi, int si = 0) const { return face_planes[si][fi]; }
183  plane_type& face_plane(int fi, int si = 0) { return face_planes[si][fi]; }
185  plane_type compute_face_plane(int fi, int si = 0) const {
186  point_type center(0, 0, 0);
187  point_type normal(0, 0, 0);
188  int last_vi = face(fi, si).back();
189  for (int vi : face(fi, si)) {
190  center += position(vi);
191  normal += cross(position(last_vi), position(vi));
192  last_vi = vi;
193  }
194  normal.normalize();
195  center *= T(1)/face(fi, si).size();
196  return plane_type(normal(0), normal(1), normal(2), -dot(center, normal));
197  }
199  int add_shell() {
200  face_planes.push_back(std::vector<plane_type>());
201  return convex_polyhedron_base::add_shell();
202  }
204  void del_shell(int si) {
205  face_planes.erase(face_planes.begin() + si);
206  convex_polyhedron_base::del_shell(si);
207  }
209  void copy_shell(int si_source, int si_target) {
210  face_planes[si_target] = face_planes[si_source];
211  convex_polyhedron_base::copy_shell(si_source, si_target);
212  }
214  void swap_shells(int si0, int si1) {
215  std::swap(face_planes[si0], face_planes[si1]);
216  convex_polyhedron_base::swap_shells(si0, si1);
217  }
219  int add_face(int si = 0) {
220  face_planes[si].push_back(plane_type());
221  return convex_polyhedron_base::add_face(si);
222  }
224  void del_face(int fi, int si = 0) {
225  face_planes[si].erase(face_planes[si].begin() + fi);
226  convex_polyhedron_base::del_face(fi, si);
227  }
229 
231  void construct_box(const box_type& box);
233  void clear();
235  void compute_vertex_locations(const plane_type& plane, std::vector<VertexPlaneLocation>& vertex_locations, T epsilon = 16 * std::numeric_limits<T>::epsilon(), std::vector<T>* vertex_signed_distances_ptr = 0) const;
237  bool check_face_locations(unsigned si, const std::vector<VertexPlaneLocation>& vertex_locations, std::vector<FacePlaneLocation>& face_locations) const;
239  bool ensure_vertex_location_consistency(unsigned si, const std::vector<T>& vertex_signed_distances, std::vector<VertexPlaneLocation>& vertex_locations, const std::vector<FacePlaneLocation>& face_locations, T epsilon = 16 * std::numeric_limits<T>::epsilon(), T epsilon_flexibiliy = 8);
241 
243  int clip_to_inside_of_plane(unsigned si, const plane_type& clip_plane, bool keep_original_shell = false, T epsilon = 16 * std::numeric_limits<T>::epsilon(), T epsilon_flexibiliy = 8);
245 
247  int clip_to_outside_of_plane(unsigned si, const plane_type& clip_plane, bool keep_original_shell = false, T epsilon = 16 * std::numeric_limits<T>::epsilon(), T epsilon_flexibiliy = 8);
249 
251  std::pair<int,int> split_at_plane(unsigned si, const plane_type& split_plane, bool keep_original_shell = false, T epsilon = 16 * std::numeric_limits<T>::epsilon(), T epsilon_flexibiliy = 8, std::vector<VertexPlaneLocation>* vertex_locations_ptr = 0);
253  std::vector<vertex_type> compute_intersection_polygon(unsigned si, const plane_type& plane, T epsilon = 16 * std::numeric_limits<T>::epsilon(), T epsilon_flexibiliy = 8) const;
254 };
255 
257 template <typename T, cgv::type::int32_type TCDim>
258 void convex_polyhedron<T,TCDim>::construct_box(const box_type& box)
259 {
260  int i, c;
261  // construct corner vertices
262  vertices.resize(8);
263  for (i = 0; i < 8; ++i) {
264  vertex_type &V = vertices[i];
265  for (c = 0; c < 3; ++c) {
266  V.position(c) = ((i & (1 << c)) == 0) ? box.get_min_pnt()(c) : box.get_max_pnt()(c);
267  if (c < texcoord_dim)
268  V.texcoord(c) = ((i & (1 << c)) == 0) ? T(0) : T(1);
269  }
270  }
271  // construct node
272  int si = add_shell();
273  shell_type& S = shell(si);
274  // construct faces
275  S.clear();
276  for (i = 0; i < 6; ++i) {
277  int fi = add_face(si);
278  plane_type& P = face_plane(fi,si);
279  P.zeros();
280  P(i / 2) = ((i & 1) == 0) ? T(1) : T(-1);
281  P(3) = ((i & 1) == 0) ? -box.get_max_pnt()(i / 2) : box.get_min_pnt()(i / 2);
282  face_type& F = face(fi, si);
283  for (unsigned j = 0; j < 4; ++j)
284  F.push_back(box_get_face_vertex_index(i, j));
285  }
286 }
287 
289 template <typename T, cgv::type::int32_type TCDim>
290 void convex_polyhedron<T, TCDim>::clear()
291 {
292  vertices.clear();
293  while (get_nr_shells() > 0)
294  del_shell((int)get_nr_shells() - 1);
295 }
296 
298 template <typename T, cgv::type::int32_type TCDim>
299 bool convex_polyhedron<T, TCDim>::check_face_locations(unsigned si, const std::vector<VertexPlaneLocation>& vertex_locations, std::vector<FacePlaneLocation>& face_locations) const
300 {
301  bool valid = true;
302  const shell_type& S = shell(si);
303  unsigned nr_faces = (unsigned)get_nr_faces(si);
304  unsigned nr_face_touches = 0;
305  face_locations.resize(nr_faces);
306  for (unsigned fi = 0; fi < nr_faces; ++fi) {
307  const face_type& F = S[fi];
308  face_locations[fi] = get_face_plane_location(face(fi, si), vertex_locations);
309  if (face_locations[fi] == FPL_TOUCH_FACE)
310  ++nr_face_touches;
311  if (!is_valid_face_location(face_locations[fi]))
312  valid = false;
313  }
314  if (nr_face_touches > 1)
315  valid = false;
316  return valid;
317 }
318 
320 template <typename T, cgv::type::int32_type TCDim>
321 bool convex_polyhedron<T, TCDim>::ensure_vertex_location_consistency(unsigned ni, const std::vector<T>& vertex_signed_distances, std::vector<VertexPlaneLocation>& vertex_locations, const std::vector<FacePlaneLocation>& face_locations, T epsilon, T epsilon_flexibiliy)
322 {
323  return false;
324 }
325 
327 template <typename T, cgv::type::int32_type TCDim>
328 int convex_polyhedron<T, TCDim>::clip_to_inside_of_plane(unsigned si, const plane_type& clip_plane, bool keep_original_shell, T epsilon, T epsilon_flexibiliy)
329 {
330  std::pair<int, int> res = split_at_plane(si, clip_plane, keep_original_shell, epsilon, epsilon_flexibiliy);
331  if (res.second == -1)
332  return res.first;
333  if (res.first == -1)
334  return -1;
335  if (res.first < res.second) {
336  del_shell(res.second);
337  return res.first;
338  }
339  copy_shell(res.first, res.second);
340  del_shell(res.first);
341  return res.second;
342 }
343 
345 template <typename T, cgv::type::int32_type TCDim>
346 int convex_polyhedron<T, TCDim>::clip_to_outside_of_plane(unsigned si, const plane_type& clip_plane, bool keep_original_shell, T epsilon, T epsilon_flexibiliy)
347 {
348  std::pair<int, int> res = split_at_plane(si, clip_plane, keep_original_shell, epsilon, epsilon_flexibiliy);
349  if (res.first == -1)
350  return res.second;
351  if (res.second == -1)
352  return -1;
353  if (res.second < res.first) {
354  del_shell(res.first);
355  return res.second;
356  }
357  copy_shell(res.second, res.first);
358  del_shell(res.second);
359  return res.first;
360 }
361 
362 
363 template <typename T, cgv::type::int32_type TCDim>
364 void convex_polyhedron<T, TCDim>::compute_vertex_locations(const plane_type& plane, std::vector<VertexPlaneLocation>& vertex_locations, T epsilon, std::vector<T>* vertex_signed_distances_ptr) const
365 {
366  // compute signed distances of vertices from clipping plane and compute sign including 0
367  unsigned vi;
368  vertex_locations.resize(get_nr_vertices());
369  std::vector<T> vertex_signed_distances_local;
370  std::vector<T>& vertex_signed_distances = vertex_signed_distances_ptr ? *vertex_signed_distances_ptr : vertex_signed_distances_local;
371  vertex_signed_distances.resize(get_nr_vertices());
372  for (vi = 0; vi < get_nr_vertices(); ++vi) {
373  T sd = dot(position(vi), reinterpret_cast<const point_type&>(plane)) + plane(3);
374  vertex_signed_distances[vi] = sd;
375  vertex_locations[vi] = (sd < -epsilon) ? VPL_INSIDE : ((sd > epsilon) ? VPL_OUTSIDE : VPL_TOUCH);
376  }
377 }
378 
380 template <typename T, cgv::type::int32_type TCDim>
381 std::pair<int, int> convex_polyhedron<T, TCDim>::split_at_plane(unsigned si, const plane_type& split_plane, bool keep_original_shell, T epsilon, T epsilon_flexibiliy, std::vector<VertexPlaneLocation>* vertex_locations_ptr)
382 {
383  // first compute vertex locations
384  std::vector<VertexPlaneLocation> vertex_locations_local;
385  std::vector<VertexPlaneLocation>& vertex_locations = vertex_locations_ptr ? *vertex_locations_ptr : vertex_locations_local;
386  std::vector<T> vertex_signed_distances;
387  compute_vertex_locations(split_plane, vertex_locations, epsilon, &vertex_signed_distances);
388 
389  // next compute face locations and check for invalid cases
390  std::vector<FacePlaneLocation> face_locations;
391  if (!check_face_locations(si, vertex_locations, face_locations)) {
392  if (!ensure_vertex_location_consistency(si, vertex_signed_distances, vertex_locations, face_locations)) {
393  std::cerr << "found case that could not be made consistent" << std::endl;
394  abort();
395  }
396  }
397 
398  // check if intersection is empty
399  VertexPlaneLocation shell_location = shell_location_side(face_locations);
400  // if complete shell is inside, no split computation is necessary
401  if (shell_location == VPL_INSIDE) {
402  // in case the original shell is not kept, reuse it as result for interior shell
403  if (!keep_original_shell)
404  return std::pair<int, int>(si, -1);
405  // otherwise construct copy of shell
406  int new_si = add_shell();
407  copy_shell(si, new_si);
408  return std::pair<int, int>(new_si, -1);
409  }
410  // if complete shell is outside
411  if (shell_location == VPL_OUTSIDE) {
412  // in case the original shell is not kept, reuse it as result for exterior shell
413  if (!keep_original_shell)
414  return std::pair<int, int>(-1, si);
415  // otherwise construct copy of shell
416  int new_si = add_shell();
417  copy_shell(si, new_si);
418  return std::pair<int, int>(-1, new_si);
419  }
420 
421  // collect halfedges of face of new plane
422  std::vector<std::pair<int, int> > new_face_halfedges;
423  // use a map to hash the newly created vertices per edge
424  std::map<std::pair<int, int>, int> edge_vertex_indices;
425  // use a map to hash touched edges with their locations to allow adding edges of the new face
426  std::map<std::pair<int, int>, FacePlaneLocation> edge_touch_events;
427 
428  // create two new shells
429  int si_inside = add_shell();
430  int si_outside = add_shell();
431  // go through all faces
432  for (unsigned fi = 0; fi < get_nr_faces(si); ++fi) {
433  switch (face_locations[fi]) {
434  // outside faces that do not touch are only touch in vertex or edge are kept
435  case FPL_OUTSIDE + FPL_TOUCH_EDGE:
436  {
437  // check touching edges for matching insides to also consider them for edges of the new face
438  std::pair<int, int> halfedge = extract_touching_halfedge(face(fi,si), vertex_locations);
439  std::pair<int, int> edge(std::min(halfedge.first, halfedge.second), std::max(halfedge.first, halfedge.second));
440  auto ete_iter = edge_touch_events.find(edge);
441  // if yes
442  if (ete_iter != edge_touch_events.end()) {
443  // check for opposite side
444  if (ete_iter->second == FPL_INSIDE + FPL_TOUCH_EDGE)
445  new_face_halfedges.push_back(std::pair<int, int>(halfedge.second, halfedge.first));
446  // erase matched and unmatched events
447  edge_touch_events.erase(ete_iter);
448  }
449  // otherwise create new edge touch event
450  else
451  edge_touch_events[edge] = FacePlaneLocation(FPL_OUTSIDE + FPL_TOUCH_EDGE);
452  }
453  case FPL_OUTSIDE + FPL_TOUCH_NOTHING:
454  case FPL_OUTSIDE + FPL_TOUCH_VERTEX:
455  {
456  int fi_outside = add_face(si_outside);
457  face(fi_outside, si_outside) = face(fi, si);
458  face_plane(fi_outside, si_outside) = face_plane(fi, si);
459  break;
460  }
461  // inside faces that do not touch are only touch in vertex or edge are kept
462  case FPL_INSIDE + FPL_TOUCH_EDGE:
463  {
464  // check touching edges for matching outsides to also consider them for edges of the new face
465  std::pair<int, int> halfedge = extract_touching_halfedge(face(fi, si), vertex_locations);
466  std::pair<int, int> edge(std::min(halfedge.first, halfedge.second), std::max(halfedge.first, halfedge.second));
467  auto ete_iter = edge_touch_events.find(edge);
468  // if yes
469  if (ete_iter != edge_touch_events.end()) {
470  // check for opposite side
471  if (ete_iter->second == FPL_OUTSIDE + FPL_TOUCH_EDGE)
472  new_face_halfedges.push_back(halfedge);
473  // erase matched and unmatched events
474  edge_touch_events.erase(ete_iter);
475  }
476  // otherwise create new edge touch event
477  else
478  edge_touch_events[edge] = FacePlaneLocation(FPL_INSIDE + FPL_TOUCH_EDGE);
479  }
480  case FPL_INSIDE + FPL_TOUCH_NOTHING:
481  case FPL_INSIDE + FPL_TOUCH_VERTEX:
482  {
483  int fi_inside = add_face(si_inside);
484  face(fi_inside, si_inside) = face(fi, si);
485  face_plane(fi_inside, si_inside) = face_plane(fi, si);
486  break;
487  }
488  // completely touching faces are kept and dublicated
489  case FPL_TOUCH_FACE:
490  {
491  int fi_outside = add_face(si_outside);
492  int fi_inside = add_face(si_inside);
493  face(fi_outside, si_outside) = face(fi_inside, si_inside) = face(fi, si);
494  face_plane(fi_outside, si_outside) = -face_plane(fi, si);
495  face_plane(fi_inside, si_inside) = face_plane(fi, si);
496  if (orientation_match(split_plane, face_plane(fi, si)))
497  change_orientation(face(fi_outside, si_outside));
498  else
499  change_orientation(face(fi_inside, si_inside));
500  break;
501  }
502  // in the split cases we need to iterate the halfedges
503  case FPL_SPLIT_TWO_VERTICES:
504  case FPL_SPLIT_VERTEX_EDGE:
505  case FPL_SPLIT_TWO_EDGES:
506  {
507  const face_type& F = face(fi, si);
508  // add an empty face on both sides
509  int fi_outside = add_face(si_outside);
510  int fi_inside = add_face(si_inside);
511  face_type& F_outside = face(fi_outside, si_outside);
512  face_type& F_inside = face(fi_inside, si_inside);
513  // copy face plane to both sides
514  face_plane(fi_outside, si_outside) = face_plane(fi_inside, si_inside) = face_plane(fi, si);
515  // loop face vertex indices and update counts
516  std::vector<int> split_vis;
517  int last_vi = F.back();
518  VertexPlaneLocation last_vertex_location = vertex_locations[last_vi];
519  for (unsigned i = 0; i < F.size(); ++i) {
520  int vi = F[i];
521  VertexPlaneLocation vertex_location = vertex_locations[vi];
522  // touching vertices go into both face loops
523  if (vertex_location == VPL_TOUCH) {
524  F_outside.push_back(vi);
525  F_inside.push_back(vi);
526  split_vis.push_back(vi);
527  if (split_vis.size() == 2 && last_vertex_location == VPL_INSIDE)
528  std::swap(split_vis[0], split_vis[1]);
529  }
530  // if previous halfedge is split, generate split vertex and add to both loops
531  else {
532  if (last_vertex_location != VPL_TOUCH && vertex_location != last_vertex_location) {
533  // check if split vertex has already been created on inverse halfedge
534  std::pair<int, int> edge(std::min(last_vi, vi), std::max(last_vi, vi));
535  auto evi_iter = edge_vertex_indices.find(edge);
536  int split_vi;
537  // if yes, get vertex index and remove entry from map
538  if (evi_iter != edge_vertex_indices.end()) {
539  split_vi = evi_iter->second;
540  edge_vertex_indices.erase(evi_iter);
541  }
542  // otherwise create new vertex on edge by affine interpolation according to unsigned distances
543  else {
544  split_vi = (int)vertices.size();
545  T distance = abs(vertex_signed_distances[vi]);
546  T last_distance = abs(vertex_signed_distances[last_vi]);
547  T lambda = last_distance / (last_distance + distance);
548  vertex_type split_vertex;
549  split_vertex.position = (T(1) - lambda) * position(last_vi) + lambda * position(vi);
550  split_vertex.texcoord = (T(1) - lambda) * texcoord(last_vi) + lambda * texcoord(vi);
551  vertices.push_back(split_vertex);
552  vertex_locations.push_back(VPL_TOUCH);
553  edge_vertex_indices[edge] = split_vi;
554  }
555  split_vis.push_back(split_vi);
556  if (split_vis.size() == 2 && last_vertex_location == VPL_INSIDE)
557  std::swap(split_vis[0], split_vis[1]);
558 
559  // add split vertex to both faces
560  F_outside.push_back(split_vi);
561  F_inside.push_back(split_vi);
562  }
563  // then add current vertex only to current loop
564  if (vertex_location == VPL_INSIDE)
565  F_inside.push_back(vi);
566  else
567  F_outside.push_back(vi);
568  }
569  last_vi = vi;
570  last_vertex_location = vertex_location;
571  }
572  if (split_vis.size() != 2) {
573  std::cerr << "expected two split vertices but found " << split_vis.size() << std::endl;
574  abort();
575  }
576  new_face_halfedges.push_back(std::pair<int, int>(split_vis[0], split_vis[1]));
577  break;
578  }
579  default:
580  std::cerr << "did not expect face plane location " << face_locations[fi] << std::endl;
581  abort();
582  }
583  }
584  // finally we need to add one face for the new plane on each side
585 
586  // prepare new faces
587  int fi_outside = add_face(si_outside);
588  int fi_inside = add_face(si_inside);
589  face_type& F_outside = face(fi_outside, si_outside);
590  face_type& F_inside = face(fi_inside, si_inside);
591  face_plane(fi_outside, si_outside) = -split_plane;
592  face_plane(fi_inside, si_inside) = split_plane;
593 
594  // first link the new face for the outside shell
595  F_outside.resize(new_face_halfedges.size());
596  F_outside[0] = new_face_halfedges[0].first;
597  int last_vi = F_outside[1] = new_face_halfedges[0].second;
598  unsigned k0 = 1;
599  for (unsigned i = 2; i < F_outside.size(); ++i) {
600  bool found = false;
601  for (unsigned k = k0; k < new_face_halfedges.size(); ++k) {
602  if (new_face_halfedges[k].first == last_vi) {
603  last_vi = F_outside[i] = new_face_halfedges[k].second;
604  if (k > k0)
605  std::swap(new_face_halfedges[k0], new_face_halfedges[k]);
606  ++k0;
607  found = true;
608  break;
609  }
610  }
611  if (!found) {
612  std::cerr << "loop of new face not complete" << std::endl;
613  abort();
614  }
615  }
616  if (k0 < new_face_halfedges.size()-1) {
617  std::cerr << "loop of new face does not use all new halfedges" << std::endl;
618  abort();
619  }
620  if (new_face_halfedges.back().first != F_outside.back() ||
621  new_face_halfedges.back().second != F_outside.front()) {
622  std::cerr << "loop of new face does is not closed by last new halfedge" << std::endl;
623  abort();
624  }
625  // finally copy to inside shell in reverse order
626  std::copy(F_outside.rbegin(), F_outside.rend(), std::back_inserter(F_inside));
627 
628  // in case original shell is kept, return pair of shell indices
629  if (keep_original_shell)
630  return std::pair<int, int>(si_inside, si_outside);
631  // otherwise overwrite original shell with outside shell and remove outside shell again
632  copy_shell(si_outside, si);
633  del_shell(si_outside);
634  return std::pair<int, int>(si_inside, si);
635 }
636 
638 template <typename T, cgv::type::int32_type TCDim>
639 std::vector<typename convex_polyhedron<T, TCDim>::vertex_type> convex_polyhedron<T, TCDim>::compute_intersection_polygon(unsigned si, const plane_type& plane, T epsilon, T epsilon_flexibiliy) const
640 {
641  // first compute vertex locations
642  std::vector<VertexPlaneLocation> vertex_locations;
643  std::vector<T> vertex_signed_distances;
644  compute_vertex_locations(plane, vertex_locations, epsilon, &vertex_signed_distances);
645 
646  // next compute face locations and check for invalid cases
647  std::vector<FacePlaneLocation> face_locations;
648  if (!check_face_locations(si, vertex_locations, face_locations)) {
649  std::cerr << "no valid face locations found" << std::endl;
650  return std::vector<vertex_type>();
651  }
652 
653  // check if intersection is empty
654  if (shell_location_side(face_locations) != VPL_TOUCH)
655  return std::vector<vertex_type>();
656 
657  // check if intersection corresponds to a touching face
658  for (unsigned fi = 0; fi < face_locations.size(); ++fi) {
659  if (face_locations[fi] == FPL_TOUCH_FACE) {
660  std::vector<vertex_type> intersection;
661  const face_type& F = face(fi, si);
662  if (orientation_match(face_plane(fi, si), plane)) {
663  for (auto v = F.begin(); v != F.end(); ++v)
664  intersection.push_back(vertices[*v]);
665  }
666  else {
667  for (auto v = F.rbegin(); v != F.rend(); ++v)
668  intersection.push_back(vertices[*v]);
669  }
670  return intersection;
671  }
672  }
673  // collect halfedges of face of new plane
674  std::vector<std::pair<int, int> > new_face_halfedges;
675  // use a map to hash the newly created vertices per edge
676  std::map<std::pair<int, int>, int> edge_vertex_indices;
677  // use a map to hash touched edges with their locations to allow adding edges of the new face
678  std::map<std::pair<int, int>, FacePlaneLocation> edge_touch_events;
679  // keep a vector of all vertices created on split edges
680  std::vector<vertex_type> split_vertices;
681 
682  // go through all faces
683  for (unsigned fi = 0; fi < get_nr_faces(si); ++fi) {
684  switch (face_locations[fi]) {
685  // outside faces can only contribute touching edges
686  case FPL_OUTSIDE + FPL_TOUCH_EDGE:
687  case FPL_INSIDE + FPL_TOUCH_EDGE:
688  {
689  // check touching edges for matching insides to also consider them for edges of the new face
690  std::pair<int, int> halfedge = extract_touching_halfedge(face(fi, si), vertex_locations);
691  std::pair<int, int> edge(std::min(halfedge.first, halfedge.second), std::max(halfedge.first, halfedge.second));
692  auto ete_iter = edge_touch_events.find(edge);
693  // if yes
694  if (ete_iter != edge_touch_events.end()) {
695  // check for opposite side
696  if (ete_iter->second != face_locations[fi]) {
697  if (ete_iter->second == FPL_OUTSIDE + FPL_TOUCH_EDGE)
698  std::swap(halfedge.first, halfedge.second);
699  new_face_halfedges.push_back(halfedge);
700  }
701  // erase matched and unmatched events
702  edge_touch_events.erase(ete_iter);
703  }
704  // otherwise create new edge touch event
705  else
706  edge_touch_events[edge] = face_locations[fi];
707  break;
708  }
709  // in the split cases we need to iterate the halfedges
710  case FPL_SPLIT_TWO_VERTICES:
711  case FPL_SPLIT_VERTEX_EDGE:
712  case FPL_SPLIT_TWO_EDGES:
713  {
714  std::vector<int> split_vis;
715  const face_type& F = face(fi, si);
716  int last_vi = F.back();
717  VertexPlaneLocation last_vertex_location = vertex_locations[last_vi];
718  for (int vi : F) {
719  VertexPlaneLocation vertex_location = vertex_locations[vi];
720  // touching vertices split the face
721  if (vertex_location == VPL_TOUCH) {
722  split_vis.push_back(vi);
723  if (split_vis.size() == 2 && last_vertex_location == VPL_OUTSIDE)
724  std::swap(split_vis[0], split_vis[1]);
725  }
726  // if previous halfedge is split, generate split vertex and add to both loops
727  else {
728  if (last_vertex_location != VPL_TOUCH && vertex_location != last_vertex_location) {
729  // check if split vertex has already been created on inverse halfedge
730  std::pair<int, int> edge(std::min(last_vi, vi), std::max(last_vi, vi));
731  auto evi_iter = edge_vertex_indices.find(edge);
732  int split_vi;
733  // if yes, get vertex index and remove entry from map
734  if (evi_iter != edge_vertex_indices.end()) {
735  split_vi = evi_iter->second;
736  edge_vertex_indices.erase(evi_iter);
737  }
738  // otherwise create new vertex on edge by affine interpolation according to unsigned distances
739  else {
740  split_vi = -1-(int)split_vertices.size();
741  T distance = abs(vertex_signed_distances[vi]);
742  T last_distance = abs(vertex_signed_distances[last_vi]);
743  T lambda = last_distance / (last_distance + distance);
744  vertex_type split_vertex;
745  split_vertex.position = (T(1) - lambda) * position(last_vi) + lambda * position(vi);
746  split_vertex.texcoord = (T(1) - lambda) * texcoord(last_vi) + lambda * texcoord(vi);
747  split_vertices.push_back(split_vertex);
748  edge_vertex_indices[edge] = split_vi;
749  }
750  split_vis.push_back(split_vi);
751  if (split_vis.size() == 2 && last_vertex_location == VPL_OUTSIDE)
752  std::swap(split_vis[0], split_vis[1]);
753  }
754  }
755  last_vi = vi;
756  last_vertex_location = vertex_location;
757  }
758  if (split_vis.size() != 2) {
759  std::cerr << "expected two split vertices but found " << split_vis.size() << std::endl;
760  abort();
761  }
762  new_face_halfedges.push_back(std::pair<int, int>(split_vis[0], split_vis[1]));
763  break;
764  }
765  case FPL_INSIDE + FPL_TOUCH_NOTHING:
766  case FPL_INSIDE + FPL_TOUCH_VERTEX:
767  case FPL_OUTSIDE + FPL_TOUCH_NOTHING:
768  case FPL_OUTSIDE + FPL_TOUCH_VERTEX:
769  break;
770  default:
771  std::cerr << "did not expect face plane location " << face_locations[fi] << std::endl;
772  abort();
773  }
774  }
775  // finally we can construct intersection polygon
776  std::vector<vertex_type> intersection;
777 
778  // first link the new face for the outside shell
779  intersection.resize(new_face_halfedges.size());
780  int fst_vi = new_face_halfedges[0].first;
781  int last_vi = new_face_halfedges[0].second;
782  intersection[0] = fst_vi < 0 ? split_vertices[-fst_vi-1] : vertices[fst_vi];
783  intersection[1] = last_vi < 0 ? split_vertices[-last_vi-1] : vertices[last_vi];
784  unsigned k0 = 1;
785  for (unsigned i = 2; i < intersection.size(); ++i) {
786  bool found = false;
787  for (unsigned k = k0; k < new_face_halfedges.size(); ++k) {
788  if (new_face_halfedges[k].first == last_vi) {
789  last_vi = new_face_halfedges[k].second;
790  intersection[i] = last_vi < 0 ? split_vertices[-last_vi-1] : vertices[last_vi];
791  if (k > k0)
792  std::swap(new_face_halfedges[k0], new_face_halfedges[k]);
793  ++k0;
794  found = true;
795  break;
796  }
797  }
798  if (!found) {
799  std::cerr << "loop of new face not complete" << std::endl;
800  abort();
801  }
802  }
803  if (k0 < new_face_halfedges.size() - 1) {
804  std::cerr << "loop of new face does not use all new halfedges" << std::endl;
805  abort();
806  }
807  if (new_face_halfedges.back().first != last_vi ||
808  new_face_halfedges.back().second != fst_vi) {
809  std::cerr << "loop of new face does is not closed by last new halfedge" << std::endl;
810  abort();
811  }
812  return intersection;
813 }
814  }
815  }
816 }
817 
818 #include <cgv/config/lib_end.h>
cgv::type::int32_type
int int32_type
this type provides an 32 bit signed integer type
Definition: standard_types.h:12
cgv::math::fvec::normalize
double normalize()
normalize the vector using the L2-Norm and return the length
Definition: fvec.h:260
cgv::media::mesh::convex_polyhedron::plane_type
cgv::math::fvec< T, 4 > plane_type
plane type
Definition: convex_polyhedron.h:136
cgv::math::fvec::zeros
void zeros()
fill the vector with zeros
Definition: fvec.h:115
cgv::media::mesh::convex_polyhedron::face_plane
const plane_type & face_plane(int fi, int si=0) const
return const reference to plane of face fi in shell si
Definition: convex_polyhedron.h:181
cgv::media::mesh::convex_polyhedron::clip_to_outside_of_plane
int clip_to_outside_of_plane(unsigned si, const plane_type &clip_plane, bool keep_original_shell=false, T epsilon=16 *std::numeric_limits< T >::epsilon(), T epsilon_flexibiliy=8)
compute the exterior of the given shell si with respect to the given clipping plane and either replac...
Definition: convex_polyhedron.h:346
cgv::math::fvec< T, 3 >
cgv::media::mesh::convex_polyhedron::vertex_type
type of a vertex
Definition: convex_polyhedron.h:141
cgv::media::mesh::convex_polyhedron::ensure_vertex_location_consistency
bool ensure_vertex_location_consistency(unsigned si, const std::vector< T > &vertex_signed_distances, std::vector< VertexPlaneLocation > &vertex_locations, const std::vector< FacePlaneLocation > &face_locations, T epsilon=16 *std::numeric_limits< T >::epsilon(), T epsilon_flexibiliy=8)
if the previous check is invalid, call this function to change the vertex signs in order to validate ...
Definition: convex_polyhedron.h:321
cgv::media::mesh::convex_polyhedron::clip_to_inside_of_plane
int clip_to_inside_of_plane(unsigned si, const plane_type &clip_plane, bool keep_original_shell=false, T epsilon=16 *std::numeric_limits< T >::epsilon(), T epsilon_flexibiliy=8)
compute the interior of the given shell si with respect to the given clipping plane and either replac...
Definition: convex_polyhedron.h:328
cgv::media::mesh::convex_polyhedron::split_at_plane
std::pair< int, int > split_at_plane(unsigned si, const plane_type &split_plane, bool keep_original_shell=false, T epsilon=16 *std::numeric_limits< T >::epsilon(), T epsilon_flexibiliy=8, std::vector< VertexPlaneLocation > *vertex_locations_ptr=0)
split the given shell into interior and exterior
Definition: convex_polyhedron.h:381
cgv::media::mesh::convex_polyhedron::del_vertex
void del_vertex(int vi)
remove a vertex
Definition: convex_polyhedron.h:158
cgv::media::mesh::convex_polyhedron::swap_shells
void swap_shells(int si0, int si1)
swap two shells
Definition: convex_polyhedron.h:214
cgv::media::mesh::convex_polyhedron::texcoord_type
cgv::math::fvec< T, TCDim > texcoord_type
texture coordinate type
Definition: convex_polyhedron.h:134
cgv::media::mesh::convex_polyhedron::add_shell
int add_shell()
add a new shell and return shell index
Definition: convex_polyhedron.h:199
cgv::media::mesh::convex_polyhedron::vertices
std::vector< vertex_type > vertices
vector to store vertices of all components in one vector
Definition: convex_polyhedron.h:147
cgv::media::mesh::convex_polyhedron::construct_box
void construct_box(const box_type &box)
construct a polyhedron representing an axis aligned box
Definition: convex_polyhedron.h:258
cgv::media::mesh::convex_polyhedron::coord_type
T coord_type
declare coordinate type
Definition: convex_polyhedron.h:128
cgv::media::mesh::convex_polyhedron::compute_face_plane
plane_type compute_face_plane(int fi, int si=0) const
compute a face plane from the vertex positions
Definition: convex_polyhedron.h:185
cgv::media::axis_aligned_box
Definition: axis_aligned_box.h:15
cgv::media::plane
Definition: plane.h:14
cgv::media::mesh::convex_polyhedron::point_type
cgv::math::fvec< T, 3 > point_type
position type
Definition: convex_polyhedron.h:132
cgv::media::mesh::convex_polyhedron
simple mesh templated with per vertex position and texture coordinates and per face plane.
Definition: convex_polyhedron.h:125
cgv::math::fvec::size
static cgv::type::uint32_type size()
return number of elements
Definition: fvec.h:149
cgv::media::axis_aligned_box::get_max_pnt
const fpnt_type & get_max_pnt() const
return a const reference to corner 7
Definition: axis_aligned_box.h:62
cgv::media::axis_aligned_box::get_min_pnt
const fpnt_type & get_min_pnt() const
return a const reference to corner 0
Definition: axis_aligned_box.h:60
cgv::media::mesh::convex_polyhedron::box_type
cgv::media::axis_aligned_box< T, 3 > box_type
type of axis aligned box
Definition: convex_polyhedron.h:138
cgv::media::mesh::convex_polyhedron::orientation_match
static bool orientation_match(const plane_type &p, const plane_type &q)
check if planes are oriented in the same normal direction
Definition: convex_polyhedron.h:177
cgv::media::mesh::convex_polyhedron::compute_intersection_polygon
std::vector< vertex_type > compute_intersection_polygon(unsigned si, const plane_type &plane, T epsilon=16 *std::numeric_limits< T >::epsilon(), T epsilon_flexibiliy=8) const
compute the intersection polygon between a shell and a plane
Definition: convex_polyhedron.h:639
cgv::media::mesh::convex_polyhedron::texcoord_dim
const cgv::type::int32_type texcoord_dim
declare dimension of texture coordinates
Definition: convex_polyhedron.h:130
cgv::media::mesh::convex_polyhedron::copy_shell
void copy_shell(int si_source, int si_target)
copy a shell onto another one
Definition: convex_polyhedron.h:209
cgv::media::mesh::convex_polyhedron::compute_vertex_locations
void compute_vertex_locations(const plane_type &plane, std::vector< VertexPlaneLocation > &vertex_locations, T epsilon=16 *std::numeric_limits< T >::epsilon(), std::vector< T > *vertex_signed_distances_ptr=0) const
construct a vector of vertex locations relative to given plane
Definition: convex_polyhedron.h:364
cgv::media::mesh::convex_polyhedron::check_face_locations
bool check_face_locations(unsigned si, const std::vector< VertexPlaneLocation > &vertex_locations, std::vector< FacePlaneLocation > &face_locations) const
compute face plane locations of the faces of one shell, store them in the vector face_locations and r...
Definition: convex_polyhedron.h:299
cgv::media::mesh::convex_polyhedron::clear
void clear()
clear data structure
Definition: convex_polyhedron.h:290
cgv
the cgv namespace
Definition: vr_calib.cxx:9
cgv::media::mesh::convex_polyhedron::del_shell
void del_shell(int si)
remove the shell of the given index, be careful as shell indices large si get implizitely decreased
Definition: convex_polyhedron.h:204
cgv::media::mesh::convex_polyhedron::face_planes
std::vector< std::vector< plane_type > > face_planes
per shell a vector to store one plane per face
Definition: convex_polyhedron.h:149
cgv::media::mesh::convex_polyhedron::face_plane
plane_type & face_plane(int fi, int si=0)
return reference to plane of face fi in shell si
Definition: convex_polyhedron.h:183