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