HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GQ_PolyDelaunay.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: Delaunay triangulation tool
7  *
8  */
9 
10 #ifndef __GQ_PolyDelaunay_h__
11 #define __GQ_PolyDelaunay_h__
12 
13 #include "GQ_API.h"
14 #include <GA/GA_Edge.h>
15 #include <UT/UT_Matrix4.h>
16 #include <UT/UT_Array.h>
17 #include <UT/UT_String.h>
18 #include <UT/UT_Vector3.h>
19 
20 #include <utility>
21 
22 class GU_Detail;
23 class GQ_Detail;
24 class GQ_Point;
25 class GQ_Edge;
26 class GQ_Face;
27 class GA_Primitive;
28 class GA_EdgeGroup;
29 class GA_PrimitiveGroup;
30 class GA_PointGroup;
31 
32 /// These internal helper classes are named generically enough that they could
33 /// potentially conflict with classes in other files. We don't want to expose
34 /// their implementation, and we can't forward declare embedded helper classes,
35 /// so we resort to wrapping them in a name space as the next best thing.
36 namespace GQ_PolyDelaunayImpl {
37  class gqTriangleSplitNode;
38  class gqFaceQueue;
39  class gqGBEdgeGroup;
40  class gqEdgeQueue;
41 }
42 
43 /// Delaunay triangulation (and refinement) of a set of 2D points.
44 ///
45 /// The algorithm used is a randomized incremental algorithm, as presented in
46 /// de Berg, van Kreveld, Overmars and Schwarzkopf. "Computational Geometry",
47 /// 2nd edition, 1999.
48 ///
49 /// It's adapted from:
50 /// Guibas, Knuth and Sharir. "Randomized incremental construction of
51 /// Delaunay and Voronoi diagrams." Algorithmica, 7:381-413, 1992.
52 ///
53 /// Most of it follows the earlier paper:
54 /// Guibas and Stolfi. "Primitives for the Manipulation of General
55 /// Subdivisions and the Computation of Voronoi Diagrams", ACM Transactions
56 /// on Graphics, 4(2):74-123, April 1985.
57 /// except for the data structure for determining which triangle a point is
58 /// in.
59 ///
60 ///
61 /// After triangulation is complete, an optional Delaunay refinement
62 /// process is applied. This follows
63 /// Shewchuk. "Delaunay Refinement Algorithms for Triangular Mesh
64 /// Generation". Computational Geometry: Theory and Applications
65 /// 22(1-3):21-74, May 2002.
66 /// http://www.cs.berkeley.edu/~jrs/papers/2dj.ps
67 ///
68 /// Shewchuk describes a number of algorithms; I've implemented the
69 /// pseudocode from pages 47-48. I used the Terminator algorithm, using
70 /// concentric circles instead of diametral lenses. I haven't implemented
71 /// one part of his algorithm: the key "terminator" innovation, his
72 /// splitPermitted routine.
73 ///
74 /// A rough overview of the algorithm can be found at:
75 /// http://www-2.cs.cmu.edu/~quake/tripaper/triangle3.html
76 ///
77 /// Shewchuk's main inspiration is this paper:
78 ///
79 /// Ruppert. "A Delaunay Refinement Algorithm for Quality 2-Dimensional
80 /// Mesh Generation", Journal of Algorithms 18(3):548-585, May 1995.
81 ///
82 ///
83 ///
84 /// TODO: this could be extended to work with 3D points under a 2D
85 /// parameterization (i.e., use texture co-ordinates instead of positions).
86 /// I'm not sure how we would handle u-wrap or v-wrap, though...
87 /// There's some discussion of this problem in
88 /// Chew, L. Paul. Guaranteed-Quality Mesh Generation for Curved
89 /// Surfaces. Proceedings of the 9th Annual Symposium on Computational
90 /// Geometry, pages 274-280. ACM, May 1993.
91 ///
92 /// Better yet, we could do what Shewchuk describes in section 7.1 of his
93 /// paper: triangulate a set of 3D primitives independently, but keeping
94 /// track of insertions on the boundaries that are shared by multiple
95 /// primitives. Apparently, it would require a way of calculating geodesic
96 /// distance between points on different primitives. I'm sure other changes
97 /// to my algorithm would be required too.
98 ///
99 ///
100 /// Notes:
101 /// * We use the GA_Edge class unconventionally here. Constraint edges
102 /// are stored as GA_Edges in a GA_EdgeGroup. The endpoint order matters:
103 /// the region to the right of p0->p1 is in the "interior" of the
104 /// constraint. The optional "primitive" member of the GA_Edge points to
105 /// the constraint primitive. However, the primitive member will be null
106 /// for constraints that have no holes (i.e., have no bridge edges). This
107 /// signals that there is no need to remove triangles from the inside of
108 /// the hole.
109 
111 {
112 public:
113  /// Some local typedefs for our helper classes to avoid having to fully
114  /// qualify the names everywhere.
115  typedef GQ_PolyDelaunayImpl::gqTriangleSplitNode gqTriangleSplitNode;
116  typedef GQ_PolyDelaunayImpl::gqFaceQueue gqFaceQueue;
117  typedef GQ_PolyDelaunayImpl::gqGBEdgeGroup gqGBEdgeGroup;
118  typedef GQ_PolyDelaunayImpl::gqEdgeQueue gqEdgeQueue;
119 
120  explicit GQ_PolyDelaunay(GU_Detail *gdp);
121  virtual ~GQ_PolyDelaunay();
122 
123  /// Set the plane where triangulation happens.
124  void setPlane(const UT_Vector3 &normal, fpreal distance);
125  /// Fit a plane to the input points. Returns true if it succeeded.
126  bool fitPlane(const GA_PointGroup *pointGroup = 0,
127  UT_Vector3 *outNormal = 0,
128  fpreal *outDistance = 0);
129 
131  { myRemoveDuplicatePoints = value; }
132 
133  /// This one is mostly here for debugging - if enabled, we don't
134  /// delete the bounding triangle that's added during construction.
136  { myKeepBoundingTriangle = value; }
137  /// If enabled, delete all edges that are outside the constrained
138  /// edges.
140  { myRemoveOutsideConstraints = value; }
141  /// Allow new points - e.g., constraint edges are split to maintain
142  /// the Delaunay property.
144  { myNewConstraintPoints = enabled; }
145  // if the new point group is set by setConstraintNewPointGroup, then
146  // setting the group name is ignored
148  { myConstraintNewPointGroupName.harden(name); }
150  { myConstraintNewPointGroup = new_pt_grp; }
151  /// Allow refinement.
153  { myRefinement = enabled; }
154  void setMaxNewPoints(int maxNewPoints)
155  { myMaxNewPoints = maxNewPoints; }
156  /// Set refinement criteria. Any triangle that fails the maximum area
157  /// test will be refined by inserting a new Steiner point.
158  /// To disable this test, set to <= 0.
159  void setMaxArea(fpreal maxArea)
160  { myMaxArea = maxArea; }
161  /// Set refinement criteria. Any triangle that fails the minimum angle
162  /// test will be refined by inserting a new Steiner point.
163  /// To disable this test, set to <= 0.
164  ///
165  /// Warning: high values can cause this to halt. Shewchuk says that
166  /// there's a proof that it works for below 21 degrees, and he finds it
167  /// works below 33, but he's had problems with higher values.
168  void setMinAngle(fpreal minAngle)
169  { myMinAngle = minAngle; }
170 
171  void triangulate(const GA_PointGroup *pointGroup,
172  const GA_PrimitiveGroup *constrainedPrimGroup);
173  void buildGeometry(bool updatePointNormals = true,
174  GA_PrimitiveGroup *out_grp = NULL);
175 
176  /// Project 3D point onto plane.
177  UT_Vector2 make2D(const UT_Vector4 &pos) const;
178  UT_Vector4 make3D(const UT_Vector2 &pos) const;
179 
180 private:
181  enum gqInsertStatus
182  {
183  GQ_Success,
184  GQ_DuplicateVertex,
185  GQ_OutsideBounding
186  };
187 
188  /// A journal of events during a point insertion. It is used to allow
189  /// the point insertion to be undone.
190  struct gqJournal
191  {
192  typedef std::pair<GA_Index,GA_Index> gqEdgePointIds;
193  typedef std::pair<gqEdgePointIds, GA_Primitive *> gqConstraint;
194 
195  /// The sequence of edge flips performed, from first to last.
196  /// The point ids are from after the flip.
197  UT_Array<gqEdgePointIds> myEdgeFlips;
198  /// The constraints that have been added and removed.
199  UT_Array<gqConstraint> myConstraintsAdded;
200  UT_Array<gqConstraint> myConstraintsRemoved;
201 
202  /// If this point was inserted on an existing edge, these are the
203  /// point ids of that original edge.
204  gqEdgePointIds myOnEdge;
205  };
206 
207  /// Set up the bounding triangle, myGqd, and the locating DAG.
208  void setup();
209  /// Set up my*Transform member variables from the plane
210  /// normal/distance.
211  void calcTransforms(const UT_Vector3 &normal, fpreal distance);
212  /// Check the constraint edges to see if any of their endpoints match
213  /// an old point, and edit them to reference the new point if they
214  /// do. The "old point" can be defined either by its position (with
215  /// *no* tolerance), or by its id.
216  ///
217  /// Very inefficient.
218  enum gq_ChangeConstraintType
219  {
220  gq_ByPos,
221  gq_ByPointId
222  };
223  void changeConstraints(gq_ChangeConstraintType type,
224  const UT_Vector2 &oldPos,
225  int oldPointId,
226  GA_Offset newPoint,
227  gqGBEdgeGroup &constrainedEdges,
228  gqJournal *journal) const;
229  /// Part of the core triangulation algorithm. Insert a new point into
230  /// the existing triangulation.
231  /// If the point is intended to lie on an existing edge (e.g.,
232  /// splitting an existing edge), the "onEdge" parameter can be used to
233  /// indicate that.
234  gqInsertStatus insertPoint(int pid, GQ_Edge *onEdge,
235  gqGBEdgeGroup &constrainedEdges,
236  gqJournal *journal);
237  /// Undo an insertion.
238  void undoInsertPoint(int pid, gqGBEdgeGroup &constrainedEdges,
239  const gqJournal &journal);
240  /// Flip a single edge. Does not update the point location data
241  /// structure.
242  void flipEdge(GQ_Edge &edge, bool updateDAG);
243  /// Part of the core triangulation algorithm. Parameter point has just
244  /// been inserted, and edge is being tested to see if it's legal
245  /// according to the Delaunay criterion; if not, the edge is flipped
246  /// and we recurse.
247  void legalizeEdge(const GQ_Point &point, GQ_Edge &edge,
248  gqGBEdgeGroup &constrainedEdges,
249  gqJournal *journal);
250  /// After triangulation, ensure that all constraints are met.
251  /// Returns false if constraints are contradictory and cannot be
252  /// enforced.
253  bool enforceConstraints(gqGBEdgeGroup &constrainedEdges);
254  /// Remove the bounding triangle
255  void removeBoundingTriangle();
256  /// Remove outside edges that aren't in the given group.
257  void removeOutsideEdges(const gqGBEdgeGroup &constrainedEdges);
258  /// Remove edges that are inside holes.
259  void removeHoles(const gqGBEdgeGroup &constrainedEdges);
260  /// Destroy all edges in queue, and recursively destroy any
261  /// neighbouring edges that are not constraints. This is Shewchuk's
262  /// "triangle eating virus" used to remove edges outside the boundary
263  /// and inside holes.
264  /// For holes, it will leave one large non-triangular face after
265  /// completion.
266  void eatEdges(gqEdgeQueue &queue,
267  const gqGBEdgeGroup &constrainedEdges);
268  /// Insert new points to satisfy refinement criteria.
269  /// Returns true if the insertion ran successfully. False otherwise.
270  bool refine(gqGBEdgeGroup &constrainedEdges);
271  /// Calculate stats related to refinement of a face.
272  void calcFaceStats(const GQ_Face &face,
273  fpreal &area,
274  const gqGBEdgeGroup *constrainedEdges,
275  fpreal *minAngle) const;
276  bool isFaceTooLarge(const GQ_Face &face) const;
277  /// Test if a given face is "bad" - i.e., excessive area, or too small
278  /// an angle.
279  bool isFaceBad(const GQ_Face &face,
280  const gqGBEdgeGroup &constrainedEdges) const;
281  /// Test if the given edge is "encroached" - i.e., if the triangles on
282  /// side have an obtuse angle opposite to the edge.
283  bool isEncroached(const GQ_Edge &edge,
284  fpreal tol = 0.05f) const;
285  /// Test if the given edge is encroached to its left side only.
286  bool isEncroachedLeft(const GQ_Edge &edge,
287  fpreal tol = 0.05f) const;
288  /// Insert a point to split an non/existing edge at the given fraction
289  /// along its length.
290  GQ_Point *insertSplitPoint(GA_Offset p0, GA_Offset p1,
291  fpreal fraction,
292  gqGBEdgeGroup &constrainedEdges,
293  GQ_Edge *onEdge,
294  bool isConstraintEdge);
295  /// Split an existing edge at somewhere-near-the-midpoint.
296  GQ_Point *splitEdge(GQ_Edge &edge,
297  gqGBEdgeGroup &constrainedEdges);
298  /// Split the given existing constraint edge using the given point.
299  /// Just updates the constraintEdges group.
300  void splitConstraint(gqGBEdgeGroup &constrainedEdges,
301  const GA_Edge &gbEdge,
302  GA_Offset newPoint,
303  gqJournal *journal) const;
304  bool splitEncroachedEdges(gqGBEdgeGroup &constrainedEdges,
305  gqEdgeQueue &edgeQueue,
306  gqFaceQueue *faceQueuePtr);
307  void updateQueuesAfterNewPoint(const GQ_Point &newPoint,
308  const gqGBEdgeGroup &constrainedEdges,
309  gqEdgeQueue &edgeQueue,
310  gqFaceQueue *faceQueue);
311  void deletePoint(int pid, gqGBEdgeGroup &constrainedEdges,
312  gqEdgeQueue *edgeQueue,
313  gqFaceQueue *faceQueue);
314 
315  /// Verify that myGqd and the locating DAG make sense.
316  void verify();
317 
318  /// Transform from 3D space to plane
319  UT_Matrix4 my3to2Transform;
320  /// Transform from plane to 3D space
321  UT_Matrix4 my2to3Transform;
322  bool myNewConstraintPoints;
323  int myMaxNewPoints;
324  int myNumNewPoints;
325  bool myRefinement;
326  bool myRemoveDuplicatePoints;
327  bool myKeepBoundingTriangle;
328  bool myRemoveOutsideConstraints;
329  fpreal myMaxArea, myMinAngle;
330  bool myHasBoundingTriangle;
331  UT_String myConstraintNewPointGroupName;
332 
333  GU_Detail *myGdp;
334  GQ_Detail *myGqd;
335  GQ_Point *myOutsideP[3];
336  GA_PointGroup *myConstraintNewPointGroup;
337  /// Root of DAG for doing point->triangle lookups
338  gqTriangleSplitNode
339  *myTriSplitRoot;
340 };
341 
342 #endif
void enableRefinement(bool enabled)
Allow refinement.
#define GQ_API
Definition: GQ_API.h:10
void setRemoveOutsideConstraints(bool value)
GLsizei const GLfloat * value
Definition: glcorearb.h:824
GQ_PolyDelaunayImpl::gqTriangleSplitNode gqTriangleSplitNode
void setMinAngle(fpreal minAngle)
void setRemoveDuplicatePoints(bool value)
GLenum GLenum GLsizei const GLuint GLboolean enabled
Definition: glcorearb.h:2539
void setMaxNewPoints(int maxNewPoints)
GA_Size GA_Offset
Definition: GA_Types.h:646
GLfloat f
Definition: glcorearb.h:1926
GQ_PolyDelaunayImpl::gqFaceQueue gqFaceQueue
GQ_PolyDelaunayImpl::gqGBEdgeGroup gqGBEdgeGroup
*get result *(waiting if necessary)*A common idiom is to fire a bunch of sub tasks at the queue
Definition: thread.h:623
GLuint const GLchar * name
Definition: glcorearb.h:786
GQ_PolyDelaunayImpl::gqEdgeQueue gqEdgeQueue
void setConstraintNewPointGroupName(const UT_String &name)
void enableNewConstraintPoints(bool enabled)
fpreal64 fpreal
Definition: SYS_Types.h:277
void setMaxArea(fpreal maxArea)
Definition: core.h:1131
void setKeepBoundingTriangle(bool value)
void setConstraintNewPointGroup(GA_PointGroup *new_pt_grp)
SIM_API const UT_StringHolder distance
type
Definition: core.h:1059
GA_API const UT_StringHolder area