Geogram  Version 1.9.1-rc
A programming library of geometric algorithms
mesh_surface_intersection_internal.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2000-2022 Inria
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *
8  * * Redistributions of source code must retain the above copyright notice,
9  * this list of conditions and the following disclaimer.
10  * * Redistributions in binary form must reproduce the above copyright notice,
11  * this list of conditions and the following disclaimer in the documentation
12  * and/or other materials provided with the distribution.
13  * * Neither the name of the ALICE Project-Team nor the names of its
14  * contributors may be used to endorse or promote products derived from this
15  * software without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27  * POSSIBILITY OF SUCH DAMAGE.
28  *
29  * Contact: Bruno Levy
30  *
31  * https://www.inria.fr/fr/bruno-levy
32  *
33  * Inria,
34  * Domaine de Voluceau,
35  * 78150 Le Chesnay - Rocquencourt
36  * FRANCE
37  *
38  */
39 
40 #ifndef GEOGRAM_MESH_MESH_SURFACE_INTERSECTION_INTERNAL
41 #define GEOGRAM_MESH_MESH_SURFACE_INTERSECTION_INTERNAL
42 
48 #include <geogram/basic/common.h>
50 #include <geogram/mesh/mesh.h>
51 #include <geogram/mesh/mesh_io.h>
52 #include <geogram/mesh/index.h>
56 
57 #include <map>
58 
59 namespace GEO {
60 
70  class MeshInTriangle : public CDTBase2d {
71  public:
72 
73  typedef exact::vec3h ExactPoint;
74 
75  /***************************************************************/
76 
84  class Edge {
85  public:
86  Edge(
87  index_t v1_in = index_t(-1),
88  index_t v2_in = index_t(-1),
89  index_t f2 = index_t(-1),
90  TriangleRegion R2 = T2_RGN_T
91  ) : v1(v1_in),
92  v2(v2_in) {
93  sym.f2 = f2;
94  sym.R2 = R2;
95  }
96  index_t v1,v2; // The two extremities of the edge
97  struct { // Symbolic information: this edge = f1 /\ f2.R2
98  index_t f2;
99  TriangleRegion R2;
100  } sym;
101  };
102 
103  /***************************************************************/
104 
112  class Vertex {
113  public:
114 
115  enum Type {
116  UNINITIALIZED, MESH_VERTEX, PRIMARY_ISECT, SECONDARY_ISECT
117  };
118 
126  geo_assert(f == M->f1_);
127  type = MESH_VERTEX;
128  mit = M;
129  init_sym(f, index_t(-1), TriangleRegion(lv), T2_RGN_T);
131  }
132 
140  MeshInTriangle* M,
141  index_t f1, index_t f2,
143  ) {
144  geo_assert(f1 == M->f1_);
145  type = PRIMARY_ISECT;
146  mit = M;
147  init_sym(f1,f2,R1,R2);
149  }
150 
157  MeshInTriangle* M, const ExactPoint& point_exact_in
158  ) {
159  type = SECONDARY_ISECT;
160  mit = M;
161  init_sym(index_t(-1), index_t(-1), T1_RGN_T, T2_RGN_T);
162  init_geometry(point_exact_in);
163  }
164 
168  Vertex() {
169  type = UNINITIALIZED;
170  mit = nullptr;
171  init_sym(index_t(-1), index_t(-1), T1_RGN_T, T2_RGN_T);
172  mesh_vertex_index = index_t(-1);
173  }
174 
179  const Mesh& mesh() const {
180  return mit->mesh();
181  }
182 
188  void print(std::ostream& out=std::cerr) const;
189 
195  std::string to_string() const {
196  std::ostringstream out;
197  print(out);
198  return out.str();
199  }
200 
201  vec2 get_UV_approx() const {
202  double u = point_exact[mit->u_].estimate();
203  double v = point_exact[mit->v_].estimate();
204  double w = point_exact.w.estimate();
205  return vec2(u/w,v/w);
206  }
207 
208  protected:
209 
216  void init_sym(
218  ) {
219  sym.f1 = f1;
220  sym.f2 = f2;
221  sym.R1 = R1;
222  sym.R2 = R2;
223  mesh_vertex_index = index_t(-1);
224  }
225 
232 
237  void init_geometry(const ExactPoint& P);
238 
239  public:
240  MeshInTriangle* mit;
241  ExactPoint point_exact; // Exact homogeneous coords using expansions
242  Type type; // MESH_VERTEX, PRIMARY_ISECT or SECONDARY_ISECT
243  index_t mesh_vertex_index; // Global mesh vertex index once created
244  struct { // Symbolic information - tri-tri isect
245  index_t f1,f2; // global facet indices in mesh
246  TriangleRegion R1,R2; // triangle regions
247  } sym;
248 #ifndef GEOGRAM_USE_EXACT_NT
249  double l; // precomputed approximated (p[u]^2 + p[v]^2) / p.w^2
250 #endif
251  };
252 
253  /***************************************************************/
254 
256 
261  const Mesh& mesh() const {
262  return mesh_;
263  }
264 
270  return exact_mesh_.target_mesh();
271  }
272 
278  void set_dry_run(bool x) {
279  dry_run_ = x;
280  }
281 
287  void save_constraints(const std::string& filename) {
288  Mesh M;
289  get_constraints(M);
290  mesh_save(M,filename);
291  }
292 
293  void begin_facet(index_t f);
294 
295  index_t add_vertex(index_t f2, TriangleRegion R1, TriangleRegion R2);
296 
297  void add_edge(
298  index_t f2,
301  );
302 
306  void commit();
307 
308 
312  void clear() override;
313 
314  protected:
318  void get_constraints(Mesh& M, bool with_edges=true) const;
319 
320  vec3 mesh_vertex(index_t v) const {
321  return vec3(mesh().vertices.point_ptr(v));
322  }
323 
324  vec3 mesh_facet_vertex(index_t f, index_t lv) const {
325  index_t v = mesh().facets.vertex(f,lv);
326  return mesh_vertex(v);
327  }
328 
329  vec2 mesh_vertex_UV(index_t v) const {
330  const double* p = mesh().vertices.point_ptr(v);
331  return vec2(p[u_], p[v_]);
332  }
333 
334  vec2 mesh_facet_vertex_UV(index_t f, index_t lv) const {
335  index_t v = mesh().facets.vertex(f,lv);
336  return mesh_vertex_UV(v);
337  }
338 
339 
340  void log_err() const {
341  std::cerr << "Houston, we got a problem (while remeshing facet "
342  << f1_ << "):" << std::endl;
343  }
344 
345  protected:
346 
347  /********************** CDTBase2d overrides ***********************/
348 
356  Sign orient2d(index_t v1,index_t v2,index_t v3) const override;
357 
369  index_t v1,index_t v2,index_t v3,index_t v4
370  ) const override;
371 
392  index_t e1, index_t i, index_t j,
393  index_t e2, index_t k, index_t l
394  ) override;
395 
402  index_t e1, index_t e2, ExactPoint& I
403  ) const;
404 
412  index_t e1, index_t e2, ExactPoint& I
413  ) const;
414 
415  public:
416  void save(const std::string& filename) const override;
417 
418  protected:
419  void begin_insert_transaction() override;
420  void commit_insert_transaction() override;
421  void rollback_insert_transaction() override;
422 
423  private:
424  MeshSurfaceIntersection& exact_mesh_;
425  const Mesh& mesh_;
426  index_t f1_;
427  index_t latest_f2_;
428  index_t latest_f2_count_;
429  coord_index_t f1_normal_axis_;
430  coord_index_t u_; // = (f1_normal_axis_ + 1)%3
431  coord_index_t v_; // = (f1_normal_axis_ + 2)%3
432  vector<Vertex> vertex_;
433  vector<Edge> edges_;
434  bool has_planar_isect_;
435  bool dry_run_;
436  mutable std::map<trindex, Sign> pred_cache_;
437  bool use_pred_cache_insert_buffer_;
438  mutable std::vector< std::pair<trindex, Sign> >
439  pred_cache_insert_buffer_;
440  };
441 
442  /*************************************************************************/
443 
450  struct IsectInfo {
451  public:
452 
457  void flip() {
458  std::swap(f1,f2);
459  A_rgn_f1 = swap_T1_T2(A_rgn_f1);
460  A_rgn_f2 = swap_T1_T2(A_rgn_f2);
461  std::swap(A_rgn_f1, A_rgn_f2);
462  B_rgn_f1 = swap_T1_T2(B_rgn_f1);
463  B_rgn_f2 = swap_T1_T2(B_rgn_f2);
464  std::swap(B_rgn_f1, B_rgn_f2);
465  }
466 
472  bool is_point() const {
473  return
474  A_rgn_f1 == B_rgn_f1 &&
475  A_rgn_f2 == B_rgn_f2 ;
476  }
477 
478  index_t f1;
479  index_t f2;
480  TriangleRegion A_rgn_f1;
481  TriangleRegion A_rgn_f2;
482  TriangleRegion B_rgn_f1;
483  TriangleRegion B_rgn_f2;
484  };
485 
486  /**********************************************************************/
487 
493  public:
494  static constexpr index_t NON_MANIFOLD = index_t(-2);
496 
508  MeshSurfaceIntersection& I, bool clear_attributes,
509  double angle_tolerance = 0.0
510  );
511 
522  void get(index_t f, index_t group_id);
523 
533 
538  void save_borders(const std::string& filename);
539 
545  void save_facet_group(const std::string& filename);
546 
552  void triangulate();
553 
554  protected:
555 
562 
571  const ExactPoint& p1, const ExactPoint& p2, const ExactPoint& p3
572  );
573 
587  const ExactPoint& P1, const ExactPoint& P2,
588  const ExactPoint& P3, const ExactPoint& P4
589  ) const;
590 
601  const ExactPoint& P1, const ExactPoint& P2, const ExactPoint& P3
602  ) const;
603 
604 
605 
606  public:
607  ExactCDT2d CDT;
608 
614  return facets_.size();
615  }
616 
623  void mark_facets(vector<index_t>& facet_is_marked) {
624  for(index_t f: facets_) {
625  facet_is_marked[f] = 1;
626  }
627  }
628 
629  private:
630  MeshSurfaceIntersection& intersection_;
631  Mesh& mesh_;
632  double angle_tolerance_;
633  index_t group_id_;
634  Attribute<index_t> facet_group_;
635  Attribute<bool> keep_vertex_;
636  Attribute<bool> c_is_coplanar_;
637  vector<bool> f_visited_;
638  vector<bool> h_visited_;
639  vector<bool> v_visited_;
640  vector<index_t> v_idx_;
641  coord_index_t u_;
642  coord_index_t v_;
643 
644  /***********************************************************/
645 
646  vector<index_t> vertices_;
647  vector<index_t> facets_;
648 
652  class Halfedges {
653  public:
654 
659  Halfedges(
660  CoplanarFacets& coplanar_facets
661  ) : mesh_(coplanar_facets.mesh_) {
662  }
663 
668  void initialize() {
669  // Use resize rather than assign so that we do not traverse
670  // all the halfedges of the mesh_
671  v_first_halfedge_.resize(mesh_.vertices.nb(), NO_INDEX);
672  h_next_around_v_.resize(mesh_.facet_corners.nb(), NO_INDEX);
673  // We only need to reset the halfedges of this set of coplanar
674  // facets.
675  for(index_t h: halfedges_) {
676  v_first_halfedge_[vertex(h,0)] = NO_INDEX;
677  h_next_around_v_[h] = NO_INDEX;
678  }
679  halfedges_.resize(0);
680  }
681 
688  index_t vertex(index_t h, index_t dlv) const {
689  index_t f = h/3;
690  index_t lv = (h+dlv)%3;
691  return mesh_.facets.vertex(f,lv);
692  }
693 
699  index_t facet(index_t h) const {
700  return h/3;
701  }
702 
709  index_t alpha2(index_t h) const {
710  index_t t1 = facet(h);
711  index_t t2 = mesh_.facet_corners.adjacent_facet(h);
712  if(t2 == NO_INDEX) {
713  return NO_INDEX;
714  }
715  for(index_t h2: mesh_.facets.corners(t2)) {
716  if(mesh_.facet_corners.adjacent_facet(h2) == t1) {
717  return h2;
718  }
719  }
721  }
722 
727  void add(index_t h) {
728  halfedges_.push_back(h);
729  index_t v1 = vertex(h,0);
730  h_next_around_v_[h] = v_first_halfedge_[v1];
731  v_first_halfedge_[v1] = h;
732  }
733 
738  vector<index_t>::const_iterator begin() const {
739  return halfedges_.begin();
740  }
741 
746  vector<index_t>::const_iterator end() const {
747  return halfedges_.end();
748  }
749 
755  index_t vertex_first_halfedge(index_t v) const {
756  return v_first_halfedge_[v];
757  }
758 
765  index_t next_around_vertex(index_t h) const {
766  return h_next_around_v_[h];
767  }
768 
774  index_t nb_halfedges_around_vertex(index_t v) const {
775  index_t result = 0;
776  for(
777  index_t h = vertex_first_halfedge(v);
778  h != NO_INDEX;
779  h = next_around_vertex(h)
780  ) {
781  ++result;
782  }
783  return result;
784  }
785 
795  index_t next_along_polyline(index_t h) const {
796  index_t v2 = vertex(h,1);
797  if(nb_halfedges_around_vertex(v2) != 1) {
798  return NO_INDEX;
799  }
800  return vertex_first_halfedge(v2);
801  }
802 
803  private:
804  Mesh& mesh_;
805  vector<index_t> halfedges_;
806  vector<index_t> v_first_halfedge_;
807  vector<index_t> h_next_around_v_;
808  } halfedges_;
809 
810 
815  class Polylines {
816  public:
817 
822  Polylines(CoplanarFacets& CF) : CF_(CF) {
823  }
824 
829  void initialize() {
830  H_.resize(0);
831  polyline_start_.resize(0);
832  polyline_start_.push_back(0);
833  }
834 
839  index_t nb() const {
840  return polyline_start_.size() - 1;
841  }
842 
847  index_as_iterator begin() const {
848  return index_as_iterator(0);
849  }
850 
856  index_as_iterator end() const {
857  return index_as_iterator(nb());
858  }
859 
865  const_index_ptr_range halfedges(index_t polyline) const {
866  geo_debug_assert(polyline < nb());
867  return const_index_ptr_range(
868  H_, polyline_start_[polyline], polyline_start_[polyline+1]
869  );
870  }
871 
875  void begin_polyline() {
876  }
877 
881  void end_polyline() {
882  polyline_start_.push_back(H_.size());
883  }
884 
891  void add_halfedge(index_t h) {
892  H_.push_back(h);
893  }
894 
900  index_t first_vertex(index_t polyline) const {
901  index_t h = H_[polyline_start_[polyline]];
902  return CF_.halfedges_.vertex(h,0);
903  }
904 
912  index_t last_vertex(index_t polyline) const {
913  index_t h = H_[polyline_start_[polyline+1]-1];
914  return CF_.halfedges_.vertex(h,1);
915  }
916 
923  index_t prev_first_vertex(index_t polyline) const {
924  if(first_vertex(polyline) != last_vertex(polyline)) {
925  return NO_INDEX;
926  }
927  index_t h = H_[polyline_start_[polyline+1]-1];
928  return CF_.halfedges_.vertex(h,0);
929  }
930 
931  private:
932  CoplanarFacets& CF_;
933  vector<index_t> H_;
934  vector<index_t> polyline_start_;
935  } polylines_;
936 
937  };
938 
939  /**********************************************************************/
940 
941 }
942 
943 #endif
Constained Delaunay triangulation in 2D.
#define geo_assert_not_reached
Sets a non reachable point in the program.
Definition: assert.h:177
#define geo_assert(x)
Verifies that a condition is met.
Definition: assert.h:149
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition: assert.h:196
Specialization of Attribute for booleans.
Definition: attributes.h:1524
Base class for constrained Delaunay triangulation.
Definition: CDT_2d.h:75
Detects and retriangulates a set of coplanar facets for MeshSurfaceIntersection.
bool edges_are_colinear(const ExactPoint &P1, const ExactPoint &P2, const ExactPoint &P3) const
Tests whether two edges are co-linear.
index_t nb_facets()
Gets the number of coplanar facets.
void mark_vertices_to_keep()
Marks the vertices that need to be kept in the simplified facets.
void save_borders(const std::string &filename)
For debugging purposes, saves border edges to a file.
bool triangles_are_coplanar(const ExactPoint &P1, const ExactPoint &P2, const ExactPoint &P3, const ExactPoint &P4) const
Tests whether two adjacent triangles are coplanar.
void save_facet_group(const std::string &filename)
For debugging purposes, saves all the facets of the group to a file.
void find_coplanar_facets()
Finds all the pairs of coplanar facets.
void triangulate()
Triangulates the kept vertices.
CoplanarFacets(MeshSurfaceIntersection &I, bool clear_attributes, double angle_tolerance=0.0)
Constructs a CoplanarFacets object associated with a MeshSurfaceIntersection.
static coord_index_t triangle_normal_axis(const ExactPoint &p1, const ExactPoint &p2, const ExactPoint &p3)
Gets the coordinate along which one can project a triangle without creating degeneracies.
void get(index_t f, index_t group_id)
Gets the set of coplanar facets from a given facet and group id.
void mark_facets(vector< index_t > &facet_is_marked)
Marks the facets.
Constrained Delaunay Triangulation with vertices that are exact points. Can be used to implement 2D C...
Definition: CDT_2d.h:1537
index_t adjacent_facet(index_t c) const
Gets the facet that a corner is adjacent to.
Definition: mesh.h:887
index_t vertex(index_t f, index_t lv) const
Gets a vertex by facet and local vertex index.
Definition: mesh.h:1054
index_range corners(index_t f) const
Gets the corners of a facet.
Definition: mesh.h:1412
void init_geometry(const ExactPoint &P)
Optimizes exact numbers in generated points and computes approximate coordinates.
const Mesh & mesh() const
Gets the mesh.
Vertex(MeshInTriangle *M, index_t f1, index_t f2, TriangleRegion R1, TriangleRegion R2)
Constructor for intersections with other facets.
void init_sym(index_t f1, index_t f2, TriangleRegion R1, TriangleRegion R2)
Initializes the symbolic information of this Vertex.
std::string to_string() const
Gets a string representation of this Vertex.
Vertex(MeshInTriangle *M, const ExactPoint &point_exact_in)
Constructor for intersections between constraints.
void print(std::ostream &out=std::cerr) const
Prints this vertex.
ExactPoint compute_geometry()
Gets the geometry of this vertex.
Vertex(MeshInTriangle *M, index_t f, index_t lv)
Constructor for macro-triangle vertices.
Meshes a single triangle with the constraints that come from the intersections with the other triangl...
void get_edge_edge_intersection(index_t e1, index_t e2, ExactPoint &I) const
Computes the intersection between two edges.
index_t create_intersection(index_t e1, index_t i, index_t j, index_t e2, index_t k, index_t l) override
Given two segments that have an intersection, create the intersection.
Sign incircle(index_t v1, index_t v2, index_t v3, index_t v4) const override
Tests the relative position of a point with respect to the circumscribed circle of a triangle.
Sign orient2d(index_t v1, index_t v2, index_t v3) const override
Tests the orientation of three vertices.
void clear() override
Mesh & target_mesh()
Gets the target mesh.
void get_edge_edge_intersection_2D(index_t e1, index_t e2, ExactPoint &I) const
Auxilliary function used by get_edge_edge_intersection() for the special case when the two edges are ...
void save_constraints(const std::string &filename)
For debugging, save constraints to a file.
const Mesh & mesh() const
Gets the readonly initial mesh.
void set_dry_run(bool x)
In dry run mode, the computed local triangulations are not inserted in the global mesh....
void save(const std::string &filename) const override
Saves this CDT to a geogram mesh file.
void commit()
Creates new vertices and new triangles in target mesh.
void get_constraints(Mesh &M, bool with_edges=true) const
For debugging, copies the constraints to a mesh.
index_t nb() const
Gets the number of (sub-)elements.
Definition: mesh.h:89
Computes surface intersections.
Mesh & target_mesh()
Gets the target mesh.
const double * point_ptr(index_t v) const
Gets a point.
Definition: mesh.h:455
Represents a mesh.
Definition: mesh.h:2693
3d vector with homogeneous coordinates
Definition: vechg.h:185
Specialization of vector for elements of type bool.
Definition: memory.h:795
Vector with aligned memory allocation.
Definition: memory.h:635
index_t size() const
Gets the number of elements.
Definition: memory.h:674
Exact predicates and constructs.
Common include file, providing basic definitions. Should be included before anything else by all head...
Classes for managing tuples of indices.
The class that represents a mesh.
Functions to load and save meshes.
Functions for computing intersections between surfacic meshes and for boolean operations.
Global Vorpaline namespace.
Definition: basic.h:55
TriangleRegion
Encodes the location of a point within a triangle.
TriangleRegion swap_T1_T2(TriangleRegion R)
Replaces T1 with T2 or T2 with T1 in a region code.
vecng< 3, Numeric::float64 > vec3
Represents points and vectors in 3d.
Definition: geometry.h:65
bool mesh_save(const Mesh &M, const std::string &filename, const MeshIOFlags &ioflags=MeshIOFlags())
Saves a mesh to a file.
geo_index_t index_t
The type for storing and manipulating indices.
Definition: numeric.h:329
Sign
Integer constants that represent the sign of a value.
Definition: numeric.h:68
void initialize(int flags=GEOGRAM_INSTALL_NONE)
Initialize Geogram.
vecng< 2, Numeric::float64 > vec2
Represents points and vectors in 2d.
Definition: geometry.h:59
geo_coord_index_t coord_index_t
The type for storing coordinate indices, and iterating on the coordinates of a point.
Definition: numeric.h:363
Stores information about a triangle-triangle intersection.
bool is_point() const
Tests whether intersection is just a point.
Symbolic computation of triangle-triangle intersection.