Geogram  Version 1.9.1
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 
451  struct IsectInfo {
452  public:
453 
458  void flip() {
459  std::swap(f1,f2);
460  A_rgn_f1 = swap_T1_T2(A_rgn_f1);
461  A_rgn_f2 = swap_T1_T2(A_rgn_f2);
462  std::swap(A_rgn_f1, A_rgn_f2);
463  B_rgn_f1 = swap_T1_T2(B_rgn_f1);
464  B_rgn_f2 = swap_T1_T2(B_rgn_f2);
465  std::swap(B_rgn_f1, B_rgn_f2);
466  }
467 
473  bool is_point() const {
474  return
475  A_rgn_f1 == B_rgn_f1 &&
476  A_rgn_f2 == B_rgn_f2 ;
477  }
478 
479  index_t f1;
480  index_t f2;
481  TriangleRegion A_rgn_f1;
482  TriangleRegion A_rgn_f2;
483  TriangleRegion B_rgn_f1;
484  TriangleRegion B_rgn_f2;
485  };
486 
487  /**********************************************************************/
488 
494  public:
495  static constexpr index_t NON_MANIFOLD = index_t(-2);
497 
509  MeshSurfaceIntersection& I, bool clear_attributes,
510  double angle_tolerance = 0.0
511  );
512 
523  void get(index_t f, index_t group_id);
524 
534 
539  void save_borders(const std::string& filename);
540 
546  void save_facet_group(const std::string& filename);
547 
553  void triangulate();
554 
555  protected:
556 
563 
572  const ExactPoint& p1, const ExactPoint& p2, const ExactPoint& p3
573  );
574 
588  const ExactPoint& P1, const ExactPoint& P2,
589  const ExactPoint& P3, const ExactPoint& P4
590  ) const;
591 
602  const ExactPoint& P1, const ExactPoint& P2, const ExactPoint& P3
603  ) const;
604 
605 
606 
607  public:
608  ExactCDT2d CDT;
609 
615  return facets_.size();
616  }
617 
624  void mark_facets(vector<index_t>& facet_is_marked) {
625  for(index_t f: facets_) {
626  facet_is_marked[f] = 1;
627  }
628  }
629 
630  private:
631  MeshSurfaceIntersection& intersection_;
632  Mesh& mesh_;
633  double angle_tolerance_;
634  index_t group_id_;
635  Attribute<index_t> facet_group_;
636  Attribute<bool> keep_vertex_;
637  Attribute<bool> c_is_coplanar_;
638  vector<bool> f_visited_;
639  vector<bool> h_visited_;
640  vector<bool> v_visited_;
641  vector<index_t> v_idx_;
642  coord_index_t u_;
643  coord_index_t v_;
644 
645  /***********************************************************/
646 
647  vector<index_t> vertices_;
648  vector<index_t> facets_;
649 
653  class Halfedges {
654  public:
655 
660  Halfedges(
661  CoplanarFacets& coplanar_facets
662  ) : mesh_(coplanar_facets.mesh_) {
663  }
664 
669  void initialize() {
670  // Use resize rather than assign so that we do not traverse
671  // all the halfedges of the mesh_
672  v_first_halfedge_.resize(mesh_.vertices.nb(), NO_INDEX);
673  h_next_around_v_.resize(mesh_.facet_corners.nb(), NO_INDEX);
674  // We only need to reset the halfedges of this set of coplanar
675  // facets.
676  for(index_t h: halfedges_) {
677  v_first_halfedge_[vertex(h,0)] = NO_INDEX;
678  h_next_around_v_[h] = NO_INDEX;
679  }
680  halfedges_.resize(0);
681  }
682 
689  index_t vertex(index_t h, index_t dlv) const {
690  index_t f = h/3;
691  index_t lv = (h+dlv)%3;
692  return mesh_.facets.vertex(f,lv);
693  }
694 
700  index_t facet(index_t h) const {
701  return h/3;
702  }
703 
710  index_t alpha2(index_t h) const {
711  index_t t1 = facet(h);
712  index_t t2 = mesh_.facet_corners.adjacent_facet(h);
713  if(t2 == NO_INDEX) {
714  return NO_INDEX;
715  }
716  for(index_t h2: mesh_.facets.corners(t2)) {
717  if(mesh_.facet_corners.adjacent_facet(h2) == t1) {
718  return h2;
719  }
720  }
722  }
723 
728  void add(index_t h) {
729  halfedges_.push_back(h);
730  index_t v1 = vertex(h,0);
731  h_next_around_v_[h] = v_first_halfedge_[v1];
732  v_first_halfedge_[v1] = h;
733  }
734 
739  vector<index_t>::const_iterator begin() const {
740  return halfedges_.begin();
741  }
742 
747  vector<index_t>::const_iterator end() const {
748  return halfedges_.end();
749  }
750 
756  index_t vertex_first_halfedge(index_t v) const {
757  return v_first_halfedge_[v];
758  }
759 
766  index_t next_around_vertex(index_t h) const {
767  return h_next_around_v_[h];
768  }
769 
775  index_t nb_halfedges_around_vertex(index_t v) const {
776  index_t result = 0;
777  for(
778  index_t h = vertex_first_halfedge(v);
779  h != NO_INDEX;
780  h = next_around_vertex(h)
781  ) {
782  ++result;
783  }
784  return result;
785  }
786 
796  index_t next_along_polyline(index_t h) const {
797  index_t v2 = vertex(h,1);
798  if(nb_halfedges_around_vertex(v2) != 1) {
799  return NO_INDEX;
800  }
801  return vertex_first_halfedge(v2);
802  }
803 
804  private:
805  Mesh& mesh_;
806  vector<index_t> halfedges_;
807  vector<index_t> v_first_halfedge_;
808  vector<index_t> h_next_around_v_;
809  } halfedges_;
810 
811 
816  class Polylines {
817  public:
818 
823  Polylines(CoplanarFacets& CF) : CF_(CF) {
824  }
825 
830  void initialize() {
831  H_.resize(0);
832  polyline_start_.resize(0);
833  polyline_start_.push_back(0);
834  }
835 
840  index_t nb() const {
841  return polyline_start_.size() - 1;
842  }
843 
848  index_as_iterator begin() const {
849  return index_as_iterator(0);
850  }
851 
857  index_as_iterator end() const {
858  return index_as_iterator(nb());
859  }
860 
866  const_index_ptr_range halfedges(index_t polyline) const {
867  geo_debug_assert(polyline < nb());
868  return const_index_ptr_range(
869  H_, polyline_start_[polyline], polyline_start_[polyline+1]
870  );
871  }
872 
876  void begin_polyline() {
877  }
878 
882  void end_polyline() {
883  polyline_start_.push_back(H_.size());
884  }
885 
892  void add_halfedge(index_t h) {
893  H_.push_back(h);
894  }
895 
901  index_t first_vertex(index_t polyline) const {
902  index_t h = H_[polyline_start_[polyline]];
903  return CF_.halfedges_.vertex(h,0);
904  }
905 
913  index_t last_vertex(index_t polyline) const {
914  index_t h = H_[polyline_start_[polyline+1]-1];
915  return CF_.halfedges_.vertex(h,1);
916  }
917 
924  index_t prev_first_vertex(index_t polyline) const {
925  if(first_vertex(polyline) != last_vertex(polyline)) {
926  return NO_INDEX;
927  }
928  index_t h = H_[polyline_start_[polyline+1]-1];
929  return CF_.halfedges_.vertex(h,0);
930  }
931 
932  private:
933  CoplanarFacets& CF_;
934  vector<index_t> H_;
935  vector<index_t> polyline_start_;
936  } polylines_;
937 
938  };
939 
940  /**********************************************************************/
941 
942 }
943 
944 #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:1420
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:2701
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.