HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_TriangleMesh.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: UT_TriangleMesh.h (UT Library, C++)
7  *
8  * COMMENTS: A class to represents a constant (unmodifiable) abstract
9  * (topological) triangulation of a 2-manifold (surface) possibly with
10  * boundary, providing the usual traversal methods. Note that no geometric
11  * information is kept in this class.
12  *
13  * Be particularly careful about what a manifold "boundary point" is! The
14  * neighbourhood of a boundary point must be homeomorphic to a half-disc
15  * (in particular, the neighbourhood minus the point itself must be connected),
16  * meaning pinched boundaries (such as what you get from two triangles sharing
17  * a single point) cannot be represented by this class.
18  *
19  */
20 
21 #ifndef __UT_TriangleMesh__
22 #define __UT_TriangleMesh__
23 
24 #include "UT_API.h"
25 #include "UT_ValArray.h"
26 #include "UT_ParallelUtil.h"
27 #include "UT_Quaternion.h"
28 #include <complex>
29 
30 // We use the usual Houdini terminology of "points", "vertices", and half-edges
31 // (hedges for short): a vertex belongs to a single triangle and is mapped to a
32 // point which exists independently. A half-edge is the share of a single
33 // triangle of an edge and is oriented consistently with triangle's winding
34 // with a source and a destination point (and vertex). A neighbouring triangle
35 // sharing the same edge has its own corresponding half-edge with the same
36 // endpoints. The two half-edges are said to be *equivalent* and point to each
37 // other through the sym() method. Boundary half-edges have negative sym().
38 // Manifold requirement dictates that each non-boundary half-edge has exactly
39 // one sym() which in turn would have that half-edge as its sym(). Furthermore,
40 // the source of a half-edge is always the destination of its sym() and
41 // vice versa.
42 //
43 // Triangles are numbered between 0 and numTriangles() - 1, points are numbered
44 // 0 to numPoints() - 1, and vertices (or equivalently half-edges) are numbered
45 // 0..numHedges() - 1. Note that numHedges() is the same as numVertices() and
46 // is 3 * numTriangles().
47 //
48 // NB. Any correspondence between elements of the mesh and the world outside
49 // should be tracked externally and should never be added to this class!
50 
52 {
53 public:
54  UT_TriangleMesh() = default;
55 
56  // Build a triangulation from a list of triangle points. The length
57  // of `tri_pts` must be a multiple of 3 and should represent a manifold
58  // surface. Edges appearing to be shared by more than two triangles are
59  // automatically treated as cut-edges. Triangle edges with the same
60  // endpoints are identified as paired half-edges, unless forced apart
61  // by passing `boundary_edge_pts`, which is an array of even length broken
62  // into pairs of consecutive point numbers each representing an edge of a
63  // triangle that should not be identified with any other triangle edges
64  // and therefore forced to be a boundary edge.
65  //
66  // The created triangle mesh will have as many triangles as indicated by
67  // tri_pts (length of tri_pts divided by 3). The "adopted" number of mesh
68  // points will be one plus the largest integer in the tri_pts array
69  // although any number in the range 0 to numPoints() - 1 that does not
70  // appear in tri_pts, has no manifestation in the actual mesh. In other
71  // words there can be unused points if tri_pts array does not cover a
72  // contiguous list starting from 0.
73  //
74  // If `reindex_points` is set to true, the points are re-indexed to to the
75  // range from 0 to number of distinct entries in `tri_pts` minus one. The
76  // index of a point will agree with its rank `tri_pts` viewed as a set:
77  // the smallest number in `tri_pts` will have rank 0, second smallest 1,
78  // etc.
79  //
80  // If `hedge_pair_class` is passed, it should have the same number of
81  // entries as `tri_pts`, each corresponding to the half-edge implicitly
82  // defined by that array in each position. Two half-edges are then allowed
83  // to be paired as sym() to each other only if they have matching
84  // hedge_pair_classes. For example, if the triangles are coming from
85  // triangulation of larger polygons, putting original hedges and those
86  // created by triangulation in different classes prevents them from being
87  // identified even when endpoints match.
88 
89  void build(const UT_IntArray &tri_pts,
90  const UT_IntArray *boundary_edge_pts = nullptr,
91  bool reindex_points = false,
92  const UT_IntArray *hedge_pair_class = nullptr);
93 
95  int numTriangles() const
96  { return int(myVertexPoint.size()) / 3; }
97 
99  int numPoints() const
100  { return int(myPointVertex.size()); }
101 
103  int numHedges() const { return int(myVertexPoint.size()); }
104 
106  int numVertices() const { return numHedges(); }
107 
108  // The jth hedge of triangle t.
110  int triangleHedge(int t, int j = 0) const
111  { return 3 * t + j % 3; }
112 
113  // The jth vertex of triangle t (same as jth hedge).
115  int triangleVertex(int t, int j = 0) const
116  { return triangleHedge(t, j); }
117 
118  // Triangle of hedge h.
120  int hedgeTriangle(int h) const { return h / 3; }
121 
122  // Triangle of vertex v.
124  int vertexTriangle(int v) const
125  { return hedgeTriangle(v); }
126 
127  // The jth point of triangle t.
129  int trianglePoint(int t, int j = 0) const
130  { return myVertexPoint(3 * t + j % 3); }
131 
133  int srcVertex(int h) const { return h; }
134 
136  int dstVertex(int h) const
137  { return (h % 3 == 2) ? h - 2 : h + 1; }
138 
140  int apxVertex(int h) const
141  { return (h % 3 == 0) ? h + 2 : h - 1; }
142 
144  int onext(int h) const { return sym(lprev(h)); }
145 
147  int vertexPoint(int v) const
148  { return myVertexPoint(v); }
149 
151  int srcPoint(int h) const
152  { return myVertexPoint(srcVertex(h)); }
153 
155  int dstPoint(int h) const
156  { return myVertexPoint(dstVertex(h)); }
157 
159  int apxPoint(int h) const
160  { return myVertexPoint(apxVertex(h)); }
161 
163  int lnext(int h) const { return dstVertex(h); }
164 
166  int lprev(int h) const { return apxVertex(h); }
167 
169  int sym(int h) const { return myHedgeSym(h); }
170 
172  int pointVertex(int p) const { return myPointVertex(p); }
173 
175  int pointHedge(int p) const { return myPointVertex(p); }
176 
178  bool isBoundaryHedge(int h) const { return sym(h) < 0; }
179 
181  bool isBoundaryPoint(int p) const
182  { return isBoundaryHedge(pointHedge(p)); }
183 
184  // Number of boundary components
186  int numBoundaries() const
187  { return int(myBoundaryRepHedges.size()); }
188 
189  // Return a first (representative) hedge for component comp.
191  int boundaryHedge(int comp) const
192  { return myBoundaryRepHedges(comp); }
193 
194  // Returns the boundary component number to which h belongs, or -1 if h
195  // is an interior hedge.
197  int hedgeBoundary(int h) const
198  { return -SYSmin(0, sym(h)) - 1; }
199 
200  // Get all the boundary component hedges in order.
201  void traceBoundary(int h0, UT_IntArray &bd) const;
202 
203  int numInteriorHedges() const
204  { return myNumInteriorHedges; }
205 
206  int numBoundaryHedges() const
207  { return numHedges() - numInteriorHedges(); }
208 
209  int numBoundaryPoints() const
210  { return numBoundaryHedges(); }
211 
212  int numInteriorPoints() const
213  { return numPoints() - numBoundaryPoints(); }
214 
215  int numEdges() const
216  { return numBoundaryHedges()
217  + numInteriorHedges() / 2; }
218 
220  { return numPoints() + numTriangles()
221  - numEdges(); }
222 
223 private:
224  UT_IntArray myVertexPoint; // Points of triangles (3 per tri)
225  UT_IntArray myHedgeSym; // Syms of triangle hedges (3 per tri)
226  UT_IntArray myPointVertex; // First hedge of each point
227 
228  // One half-edge representing each boundary component.
229  UT_IntArray myBoundaryRepHedges;
230  int myNumInteriorHedges = 0;
231 };
232 
233 // An embedded triangle mesh will support a point position for each point
234 // in the mesh. The entire implementation is left to the template parameter.
235 // This includes the return type of the point position.
236 
237 template <typename E>
239 {
240 public:
241 
242  using PointCoords = typename E::value_type;
243  using Real = typename PointCoords::value_type;
244  using Complex = typename std::complex<Real>;
246 
247  explicit UT_EmbeddedTriangleMesh() = default;
248 
250  auto pointPosition(int pt) const
251  { return myEmbedding->pointPosition(pt); }
252 
253  void setMeshPointPositions(const E *ptpos)
254  { myEmbedding = ptpos; }
255 
256  bool hasEmbedding() const
257  { return myEmbedding != nullptr; }
258 
259  // Returns the barycenter of t.
261  auto triangleBarycenter(int t) const;
262 
263  // Returns the normal of the triangle t. If normalize is set to false,
264  // the normal vector will have magnitude equal to twice the area of t.
265  //
266  // NB. The normals used here are based on COUNTERCLOCKWISE orientation
267  // of triangles UNLIKE Houdini polygon normals which are backwards. Based
268  // on what these meshes are intended for, it's sheer asininity to adopt
269  // Houdini's convention here and people in disagreement should take/retake
270  // a bunch of math classes.
271 
273  auto triangleNormal(int t, bool normalize = true) const;
274 
275  // Returns the area of triangle t.
277  auto triangleArea(int t) const
278  { return 0.5 * triangleNormal(t, false).length(); }
279 
280  // Returns the (integral of) gradient of the linear function defined on a
281  // triangle t by values f0, f1, and f2 on its points 0, 1, and 2,
282  // respectively, over t.
283  template <typename T>
284  auto triangleGradient(int t, T f0, T f1, T f2) const;
285 
286  // Returns the (integral of) gradient of the function defined on dual
287  // vertices (triangle centroids) over the dual cell of a point p.
288  template <typename F, typename A>
289  Vector3 pointGradient(int p,
290  const F &tri_fn,
291  const A &pt_atlas) const;
292 
293  // Calculates the (integral of) Laplacian of a function defined on dual
294  // vertices over triangle t. The function value at triangles is supplied by
295  // tri_fn.
296  template <typename F>
297  auto triangleLaplacian(int t, const F &tri_fn) const;
298 
299  // Calculates the (integral of) Laplacian of a function defined on dual
300  // vertices over the dual cell of point p. The function value at a point
301  // is supplied by pt_fn.
302  template <typename F>
303  auto pointLaplacian(int p, const F &pt_fn) const;
304 
305  // Returns the vector from the src to dst of half-edge h.
307  auto hedgeVector(int h) const;
308 
309  // Returns the length of a half-edge h.
311  auto hedgeLength(int h) const
312  { return hedgeVector(h).length(); }
313 
314  // The cotan-based estimated length of the dual to h.
316  auto dualHedgeLength(int h) const;
317 
318  // Returns the dihedral angle between the two triangles incident to h.
319  // This is angles between the normals with a positive sign if the edge
320  // convex and negative otherwise.
321  //
322  // NB. These dihedral angles are calculating based on counterclockwise
323  // triangle normals. So, the sign must be flipped if the mesh is built
324  // consistently with a Houdini polygon geometry to give the expected
325  // dihedral angle. (See notes above `triangleNormal()`).
327  auto hedgeDihedralAngle(int h) const;
328 
329  // Returns the cotan of the opposite angle to h.
331  auto hedgeCotan(int h) const;
332 
333  // Returns the angle (in radians) of the angle opposite to h.
335  auto hedgeAngle(int h) const
336  { return SYSatan2(Real(1.0), hedgeCotan(h)); }
337 
338  // Returns the dual area of a point which is the a third of the sum of
339  // areas of triangles incident to it.
341  auto pointDualArea(int p) const;
342 
343  // Returns the sum of angles of vertices incident to p.
345  auto pointAngleSum(int p) const;
346 
347  // Returns a unit normal for a point p. By default the point normal is the
348  // normalized angle-weighted sum of the normals of the incident triangles.
349  // If `area_weighted` is set to true, incident triangle areas are used
350  // instead of angles.
351  //
352  // NB. Like triangle normals, point normals are based on counterclockwise
353  // winding of triangles. See the notes above `triangleNormal()`.
355  auto pointNormal(int p, bool area_weighted = false) const;
356 
357  // Returns the integral of the gaussian curvature in the mesh dual cell
358  // of a point p.
360  auto pointGaussCurvature(int p) const;
361 
362  // Returns the integral of the mean curvature in the mesh dual cell
363  // of a point p.
365  auto pointMeanCurvature(int p) const;
366 
367  // Returns the integral of the gaussian curvature over a mesh triangle.
368  // Note that the curvature is calculated based on the Gaussian curvature
369  // of the incident vertices and therefore depends on all the angles
370  // incident to three corners of the tringle.
372  auto triangleGaussCurvature(int t) const;
373 
374  // Returns the integral of the mean curvature in the mesh dual cell
375  // of a point p.
377  auto triangleMeanCurvature(int t) const;
378 
379  // Returns maximum principal direction as a complex number. The complex
380  // plane based on which the number should be interpreted is the plane
381  // of the triangle with the triangles first half-edge as the real axis.
384 
385 
386  // Given the polar coordinates of a tangent vector in the exponential
387  // mapping of the star of point p, generate an extrinsic (vector 3)
388  // tangent vector to the embedding of the mesh.
390  Real angular) const;
391 
392  // Same as the abvoe method, using a complex number.
395 
396 private:
397  const E *myEmbedding = nullptr;
398 };
399 
400 template <typename E>
401 auto
403 {
404  return (pointPosition(trianglePoint(t, 0))
405  + pointPosition(trianglePoint(t, 1))
406  + pointPosition(trianglePoint(t, 2))) / 3.0;
407 }
408 
409 template <typename E>
410 auto
412 {
413  auto a = pointPosition(trianglePoint(t, 0));
414  auto b = pointPosition(trianglePoint(t, 1));
415  auto c = pointPosition(trianglePoint(t, 2));
416 
417  auto nml = cross(b - a, c - a);
418  if (normalize)
419  nml.normalize();
420 
421  return nml;
422 }
423 
424 template <typename E>
425 template <typename T>
426 auto
428 {
429  auto a = pointPosition(trianglePoint(t, 0));
430  auto b = pointPosition(trianglePoint(t, 1));
431  auto c = pointPosition(trianglePoint(t, 2));
432 
433  auto nml = cross(b - a, c - a);
434  auto areax2 = nml.normalize();
435 
436  return (f2 * cross(nml, b - a) + f0 * cross(nml, c - b)
437  + f1 * cross(nml, a - c)) / areax2;
438 }
439 
440 template <typename E>
441 template <typename F>
442 auto
444 {
445  auto f = tri_fn(t);
446  auto res = 0.0 * f;
447 
448  for (int i = 0; i < 3; i++)
449  {
450  auto h = triangleHedge(t, i);
451  auto h_sym = sym(h);
452  if (h_sym < 0)
453  continue;
454 
455  auto t_sym = hedgeTriangle(h_sym);
456  auto f_sym = tri_fn(t_sym);
457 
458 #if 0
459  auto wt = 1.0 / (hedgeCotan(h) + hedgeCotan(h_sym));
460  res += (f_sym - f) * wt;
461 #else
462  res += (f_sym - f);
463 #endif
464 
465  }
466  return 0.5 * res;
467 }
468 
469 template <typename E>
470 template <typename F>
471 auto
473 {
474  auto f = pt_fn(p);
475 
476  auto res = 0.0 * f;
477  auto ring_area = 0.0;
478 
479  auto h0 = pointHedge(p);
480  auto h = h0;
481  auto h_sym = -1;
482  do
483  {
484  // auto cot = SYSmax(1e-7, hedgeCotan(h));
485  auto cot = hedgeCotan(h);
486  h_sym = sym(h);
487  if (h_sym >= 0)
488  // cot += SYSmax(1e-7, hedgeCotan(h_sym));
489  cot += hedgeCotan(h_sym);
490 
491 
492  auto f_dst = pt_fn(dstPoint(h));
493  res += cot * (f_dst - f);
494  ring_area += triangleArea(hedgeTriangle(h));
495  } while (h = onext(h), h >= 0 && h != h0);
496 
497  if (h < 0)
498  {
499  h_sym = lprev(sym(h_sym));
500  auto cot = hedgeCotan(h_sym);
501  auto f_dst = pt_fn(srcPoint(h_sym));
502  res += cot * (f_dst - f);
503  }
504 
505  return 0.5 * res;
506 }
507 
508 template <typename E>
509 template <typename F, typename A>
512  const A &pt_atlas) const
513 {
514  auto hedge_contribution = [&](int h)
515  {
516  auto i = Complex(0.0, 1.0);
517  auto t = hedgeTriangle(h);
518  auto s = tri_fn(t);
519  auto c = hedgeCotan(h);
520  auto h_z = pt_atlas.hedgeComplexCoord(h);
521  auto z = s * c * h_z * i;
522 
523  auto h_sym = sym(h);
524  if (h_sym >= 0)
525  {
526  auto t_sym = hedgeTriangle(h_sym);
527  auto s_sym = tri_fn(t_sym);
528  z -= s_sym * c * h_z * i;
529  }
530 
531  return z;
532  };
533 
534  Complex z(0, 0);
535  auto h0 = pointHedge(p);
536  auto h = h0;
537  do
538  {
539  z += hedge_contribution(h);
540  } while (h = onext(h), h >= 0 && h != h0);
541 
542  return pt_atlas.extrinsicPointTangent(p, z) / pointDualArea(p);
543 }
544 
545 template <typename E>
546 auto
548 {
549  return (pointPosition(dstPoint(h)) - pointPosition(srcPoint(h)));
550 }
551 
552 template <typename E>
553 auto
555 {
556  auto cotsum = hedgeCotan(h);
557  auto h_sym = sym(h);
558  if (h_sym >= 0)
559  cotsum += hedgeCotan(h_sym);
560 
561  return 0.5 * hedgeLength(h) * cotsum;
562 }
563 
564 template <typename E>
565 auto
567 {
568  auto h_sym = sym(h);
569  if (h_sym < 0)
570  return Real(0.0);
571 
572  auto n0 = triangleNormal(hedgeTriangle(h));
573  auto n1 = triangleNormal(hedgeTriangle(h_sym));
574  auto e = hedgeVector(h);
575  e.normalize();
576  return SYSatan(dot(e, cross(n0, n1)), dot(n0, n1));
577 }
578 
579 template <typename E>
580 auto
582 {
583  auto a = pointPosition(dstPoint(lnext(h)));
584  auto b = pointPosition(srcPoint(h));
585  auto c = pointPosition(dstPoint(h));
586 
587  auto u = b - a;
588  auto v = c - a;
589 
590  return dot(u, v) / cross(u, v).length();
591 }
592 
593 template <typename E>
594 auto
596 {
597  auto h = pointHedge(p);
598  auto h0 = h;
599  auto sum = triangleArea(hedgeTriangle(h));
600 
601  while (h = onext(h), h >= 0 && h != h0)
602  sum += triangleArea(hedgeTriangle(h));
603 
604  return sum / Real(3.0);
605 }
606 
607 template <typename E>
608 auto
610 {
611  auto h = pointHedge(p);
612  auto h0 = h;
613  auto sum = hedgeAngle(lnext(h));
614 
615  while (h = onext(h), h >= 0 && h != h0)
616  sum += hedgeAngle(lnext(h));
617 
618  return sum;
619 }
620 
621 template <typename E>
622 auto
623 UT_EmbeddedTriangleMesh<E>::pointNormal(int p, bool area_weighted) const
624 {
625  auto h = pointHedge(p);
626  auto h0 = h;
627 
628  if (area_weighted)
629  {
630  auto nml = triangleNormal(hedgeTriangle(h), false);
631  while (h = onext(h), h >= 0 && h != h0)
632  nml += triangleNormal(hedgeTriangle(h), false);
633 
634  nml.normalize();
635  return nml;
636  }
637 
638  auto nml = triangleNormal(hedgeTriangle(h)) * hedgeAngle(lnext(h));
639  while (h = onext(h), h >= 0 && h != h0)
640  nml += triangleNormal(hedgeTriangle(h)) * hedgeAngle(lnext(h));
641 
642  nml.normalize();
643  return nml;
644 }
645 
646 template <typename E>
647 auto
649 {
650  return (isBoundaryPoint(p) ? Real(M_PI) : Real(2.0 * M_PI))
651  - pointAngleSum(p);
652 }
653 
654 template <typename E>
655 auto
657 {
658  auto h0 = pointHedge(p);
659  auto h = h0;
660  auto b = Real(0.0);
661  do
662  {
663  b += hedgeLength(h) * hedgeDihedralAngle(h);
664  } while (h = onext(h), h >= 0 && h != h0);
665  return Real(0.25) * b;
666 }
667 
668 template <typename E>
669 auto
671 {
672  auto K = Real(-M_PI);
673  for (int i = 0; i < 3; i++)
674  {
675  int h = triangleHedge(t, i);
676  int p = apxPoint(h);
677  auto a = hedgeAngle(h);
678  K += (isBoundaryPoint(p) ? 1 : 2) * M_PI * a / pointAngleSum(p);
679  }
680 
681  return K;
682 }
683 
684 template <typename E>
685 auto
687 {
688  auto H = Real(0);
689  for (int i = 0; i < 3; i++)
690  {
691  auto h = triangleHedge(t, i);
692  H += 0.25 * hedgeLength(h) * hedgeDihedralAngle(h);
693  }
694  return H;
695 }
696 
697 template<typename E>
700 {
701  Complex z(0, 0);
702  auto partial_angle_sum = Real(0.0);
703 
704  for (int i = 0; i < 3; i++)
705  {
706  auto h = triangleHedge(t, i);
707  auto l = hedgeLength(h);
708  auto alpha = hedgeDihedralAngle(h);
709 
710  auto theta = partial_angle_sum;
711 
712  // Multiply by theta by 2.0 to go from vectors to directions.
713  Complex r(SYScos(2.0 * theta), SYSsin(2.0 * theta));
714 
715  z += l * alpha * r;
716  partial_angle_sum += (M_PI - hedgeAngle(lprev(h)));
717  }
718 
719  return std::sqrt(z);
720 }
721 
722 template<typename M>
725 {
726  return extrinsicTriangleTangent(p, std::abs(z), std::arg(z));
727 }
728 
729 template<typename M>
732  Real angular) const
733 {
734  auto n = triangleNormal(t);
735  auto x = hedgeVector(triangleHedge(t, 0));
736  x.normalize();
737  x *= radial;
738  auto Jx = cross(n, x);
739  return x * SYScos(angular) + Jx * SYSsin(angular);
740 }
741 
742 
743 // UT_PointAtlas is constructs a local polar coordinate system around each
744 // point of an embedded triangle mesh to support computations on a complex
745 // plane corresponding to (the exponential map of) the ring of each point
746 // in which the real axis correspond to the half-edge returned by
747 // `pointHedge()` method of the mesh.
748 
749 template <typename M>
751 {
752 public:
753  using Real = typename M::Real;
754  using Complex = typename M::Complex;
755  using Vector3 = typename M::Vector3;
756 
757  explicit UT_PointAtlas(const M &mesh);
758 
759 
760  // Returns the exponential map scale factor at p which is 2π divided by the
761  // sum of the mesh angles incident to p.
763  Real pointAngleScale(int p) const
764  { return myPointScale(p); }
765 
766  // Returns the *raw* angle sum between h0 = mesh.pointHedge(p), where p
767  // is the source point of h, and h in a counterclockwise scale (using
768  // mesh.onext()).
770  Real hedgeAngularOffset(int h) const;
771 
772  // Returns the angular component of the polar coordinates of the vector
773  // represented by h in the exponential mapping of star of the source point
774  // of h onto a planar patch. This is the angularOffset of h scaled by point
775  // angle scale at the source of h. The output will be in range [0, 2π).
778  { return pointAngleScale(myMesh.srcPoint(h))
779  * hedgeAngularOffset(h); }
780 
781  // Returns the radial component of the polar coordinates of the vector
782  // represented by h (or equivalently the destination point of h) in the
783  // exponential mapping of star of the source point of h onto a planar
784  // patch. This is simply the length of h!
786  Real radialCoord(int h) const
787  { return myMesh.hedgeLength(h); }
788 
789  // Returns the polar coordinates of the vector represented by h in the
790  // exponential mapping of the star of the source point of h as a complex
791  // number with the same polar coordinates.
793  Complex hedgeComplexCoord(int h) const;
794 
796  Complex hedgeComplexDirection(int h, int n = 1) const;
797 
798  // Returns the the angular polar component of the complex number z. This
799  // is the phase of z, implemented using std::arg(z), in range [0, 2π).
802 
803  // Returns the radial component of the complex number z.
806  { return std::abs(z); }
807 
808  // Returns the outgoing half-edge at p with the largest angular coordinate
809  // less than or equal to `angle`.
810  //
811  // Note: This is done using a linear scan. We could be brought down to
812  // O(log(valence(p))) if we maintained a random-access list of incident
813  // hedge angular values. Since we like to be able to lookup hedge angular
814  // coordinates in constant time, this would require de facto repeated
815  // storage of the values with a rearrangement. Given that average point
816  // valence in a triangle mesh is 6, this would be excessive.
817  int angularLowerBoundHedge(int p, Real angle) const;
818 
819  // Given the polar coordinates of a tangent vector in the exponential
820  // mapping of the star of point p, generate an extrinsic (vector 3)
821  // tangent vector to the embedding of the mesh.
822  Vector3 extrinsicPointTangent(int p, Real radial,
823  Real angular) const;
824 
825  // Same as the abvoe method, using a complex number.
828 
831 
832 private:
833  const M &myMesh;
834  UT_Array<Real> myPointScale;
835  UT_Array<Real> myHedgeAngle;
836 };
837 
838 template<typename M>
840  myMesh(mesh)
841 {
842  myPointScale.setSizeNoInit(mesh.numPoints());
843  myHedgeAngle.setSizeNoInit(mesh.numHedges());
844 
845  UTparallelFor(UT_BlockedRange<int>(0, mesh.numPoints()),
846  [&](const UT_BlockedRange<int> &range)
847  {
848  for (int p = range.begin(), pe = range.end(); p != pe; p++)
849  {
850  auto h0 = mesh.pointHedge(p);
851  auto h = h0;
852  Real angle_sum = Real(0);
853  do
854  {
855  // Note that myHedgeAngle(h) carries the sum of source angles
856  // from h0 to h (inclusive). The last of these should be
857  // excluded when calculating the angle coordinate of h.
858  // We include it here in the sum because this is the only
859  // place to keep track that angle semantically (think of
860  // boundary points).
861  angle_sum += mesh.hedgeAngle(mesh.lnext(h));
862  myHedgeAngle(h) = angle_sum;
863  } while (h = mesh.onext(h), h >= 0 && h != h0);
864  myPointScale(p) = 2 * M_PI / angle_sum;
865  }
866  });
867 }
868 
869 template<typename M>
870 typename UT_PointAtlas<M>::Real
872 {
873  if (h == myMesh.pointHedge(myMesh.srcPoint(h)))
874  return 0.0;
875 
876  return myHedgeAngle(myMesh.lnext(myMesh.sym(h)));
877 }
878 
879 template<typename M>
882 {
883  auto a = hedgeAngularCoord(h);
884  return radialCoord(h) * Complex(SYScos(a), SYSsin(a));
885 }
886 
887 template<typename M>
890 {
891  auto a = hedgeAngularCoord(h);
892  return Complex(SYScos(Real(n) * a), SYSsin(Real(n) * a));
893 }
894 
895 template<typename M>
896 typename UT_PointAtlas<M>::Real
898 {
899  auto t = std::arg(z);
900  return t - 2.0 * M_PI * SYSfloor(t / (2.0 * M_PI));
901 }
902 
903 template<typename M>
904 int
906 {
907  // This can in theory be made to run in time logarithmic in valence
908  // of p; but this requires random access into the orderd list incident
909  // outgoing hedges at p. Unless we want to do lots of these queries on
910  // very irregular meshes this would be excessive.
911 
912  angle = angle - 2.0 * M_PI * SYSfloor(angle / (2.0 * M_PI));
913  angle /= pointAngleScale(p);
914 
915  auto h0 = myMesh.pointHedge(p);
916  auto h = h0;
917  auto last_h = h0;
918  do
919  {
920  if (myHedgeAngle(h) > angle)
921  return h;
922  last_h = h;
923  } while (h = myMesh.onext(h), h >= 0 && h != h0);
924 
925  // We have zoned the angle above to [0..2π) range, therefore getting
926  // here means that when dividing by pointAngleScale(p) our corrected
927  // angle ended up going over the angle sum at p which could only be
928  // due to numerical error. So, we just drop down to the last outgoing
929  // half-edge before getting back to h0.
930  return last_h;
931 }
932 
933 template<typename M>
936 {
937  auto h = angularLowerBoundHedge(p, angular);
938  auto n = myMesh.triangleNormal(myMesh.hedgeTriangle(h));
939  auto a = (angular - hedgeAngularCoord(h)) / pointAngleScale(p);
940 
941  auto v = myMesh.hedgeVector(h);
942  v.normalize();
943  v *= radial;
944 
945  if (a > 0.0)
946  {
948  q.updateFromAngleAxis(a, n);
949  v = q.rotate(v);
950  }
951  return v;
952 }
953 
954 template<typename M>
957 {
958  return extrinsicPointTangent(p, complexToRadial(z), complexToAngular(z));
959 }
960 
961 template<typename M>
964 {
965  Complex z(0, 0);
966  auto h0 = myMesh.pointHedge(p);
967  auto h = h0;
968  do
969  {
970  auto l = myMesh.hedgeLength(h);
971  auto alpha = myMesh.hedgeDihedralAngle(h);
972  z += l * alpha * hedgeComplexDirection(h, 2);
973  } while (h = myMesh.onext(h), h >= 0 && h != h0);
974 
975  return std::sqrt(z);
976 }
977 
978 #endif
979 
GLdouble s
Definition: glew.h:1390
typename M::Real Real
SYS_FORCE_INLINE int lnext(int h) const
SYS_FORCE_INLINE int apxVertex(int h) const
SYS_FORCE_INLINE int boundaryHedge(int comp) const
SYS_FORCE_INLINE int vertexPoint(int v) const
GLenum GLint * range
Definition: glew.h:3500
SYS_FORCE_INLINE GA_Offset srcPoint(const GA_Detail *gdp, GEO_Hedge h)
Definition: GEO_Hedge.h:186
SYS_FORCE_INLINE auto pointNormal(int p, bool area_weighted=false) const
GLclampf GLclampf GLclampf alpha
Definition: glew.h:1520
UT_PointAtlas(const M &mesh)
auto triangleLaplacian(int t, const F &tri_fn) const
SYS_FORCE_INLINE auto hedgeDihedralAngle(int h) const
SYS_FORCE_INLINE int triangleVertex(int t, int j=0) const
SYS_FORCE_INLINE int pointVertex(int p) const
typename M::Vector3 Vector3
SYS_FORCE_INLINE Real hedgeAngularOffset(int h) const
SYS_FORCE_INLINE Real pointAngleScale(int p) const
SYS_FORCE_INLINE auto pointAngleSum(int p) const
SYS_FORCE_INLINE Complex hedgeComplexDirection(int h, int n=1) const
SYS_FORCE_INLINE int srcVertex(int h) const
void setSizeNoInit(exint newsize)
Definition: UT_Array.h:498
internal::named_arg< T, char > arg(string_view name, const T &arg)
Definition: core.h:1393
Vector3 pointGradient(int p, const F &tri_fn, const A &pt_atlas) const
GLboolean GLboolean GLboolean GLboolean a
Definition: glew.h:9477
vfloat4 sqrt(const vfloat4 &a)
Definition: simd.h:7231
int numInteriorPoints() const
GA_Offset srcVertex(GEO_Hedge)
Definition: GEO_Hedge.h:179
#define UT_API
Definition: UT_API.h:13
SYS_FORCE_INLINE int numTriangles() const
GLdouble l
Definition: glew.h:9122
const GLdouble * v
Definition: glew.h:1391
3D Vector class.
GLdouble angle
Definition: glew.h:9135
uint64 value_type
Definition: GA_PrimCompat.h:29
#define M_PI
Definition: ImathPlatform.h:51
SYS_FORCE_INLINE auto hedgeLength(int h) const
UT_Vector3T< Real > Vector3
int numBoundaryPoints() const
SYS_FORCE_INLINE int hedgeTriangle(int h) const
GLdouble GLdouble z
Definition: glew.h:1559
Vector3 extrinsicPointTangent(int p, Real radial, Real angular) const
SYS_FORCE_INLINE int numPoints() const
SYS_FORCE_INLINE Real radialCoord(int h) const
GLdouble GLdouble GLdouble GLdouble q
Definition: glew.h:1414
void setMeshPointPositions(const E *ptpos)
SYS_FORCE_INLINE int apxPoint(int h) const
auto pointLaplacian(int p, const F &pt_fn) const
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER T abs(T a)
Definition: ImathFun.h:55
void updateFromAngleAxis(T angle, const UT_Vector3T< T > &axis, int normalize=1)
GLclampf f
Definition: glew.h:3499
GLint GLint GLint GLint GLint x
Definition: glew.h:1252
SIM_DerScalar length() const
SYS_FORCE_INLINE int numVertices() const
SYS_FORCE_INLINE int numBoundaries() const
int angularLowerBoundHedge(int p, Real angle) const
SYS_FORCE_INLINE Complex hedgeComplexCoord(int h) const
SYS_FORCE_INLINE auto pointGaussCurvature(int p) const
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:218
#define SYS_FORCE_INLINE
Definition: SYS_Inline.h:45
SYS_FORCE_INLINE Complex pointMaxPrincipalDirection(int p) const
UT_EmbeddedTriangleMesh()=default
Vector3 extrinsicTriangleTangent(int p, Real radial, Real angular) const
GLsizei n
Definition: glew.h:4040
const GLfloat * c
Definition: glew.h:16296
SYS_FORCE_INLINE bool isBoundaryPoint(int p) const
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
SYS_FORCE_INLINE Real complexToRadial(Complex z) const
SYS_FORCE_INLINE int sym(int h) const
GLfloat GLfloat GLfloat GLfloat h
Definition: glew.h:8011
SYS_API fpreal32 SYSfloor(fpreal32 val)
SYS_FORCE_INLINE int lprev(int h) const
GridType::Ptr normalize(const GridType &grid, bool threaded, InterruptT *interrupt)
Normalize the vectors of the given vector-valued grid.
SYS_FORCE_INLINE auto hedgeAngle(int h) const
typename GU_TriangleMeshDetailLinkT< T >::value_type PointCoords
int numBoundaryHedges() const
SYS_FORCE_INLINE GA_Offset dstVertex(T &iface, GEO_Hedge h)
Definition: GEO_Hedge.h:215
SYS_FORCE_INLINE Real hedgeAngularCoord(int h) const
SYS_FORCE_INLINE auto triangleArea(int t) const
SYS_FORCE_INLINE int onext(int h) const
int eulerCharacteristic() const
auto triangleGradient(int t, T f0, T f1, T f2) const
GLdouble GLdouble GLdouble b
Definition: glew.h:9122
GLfloat GLfloat p
Definition: glew.h:16321
SYS_FORCE_INLINE auto triangleMeanCurvature(int t) const
Quaternion class.
Definition: UT_Quaternion.h:51
SYS_FORCE_INLINE auto pointPosition(int pt) const
SYS_FORCE_INLINE int dstPoint(int h) const
SYS_FORCE_INLINE Complex triangleMaxPrincipalDirection(int t) const
SYS_FORCE_INLINE int numHedges() const
SYS_FORCE_INLINE auto pointMeanCurvature(int p) const
SYS_FORCE_INLINE int triangleHedge(int t, int j=0) const
SYS_FORCE_INLINE bool isBoundaryHedge(int h) const
GLdouble GLdouble GLdouble r
Definition: glew.h:1406
SYS_FORCE_INLINE auto triangleGaussCurvature(int t) const
SYS_FORCE_INLINE int vertexTriangle(int v) const
int numEdges() const
SYS_FORCE_INLINE int pointHedge(int p) const
typename M::Complex Complex
SYS_FORCE_INLINE Real complexToAngular(Complex z) const
SYS_FORCE_INLINE auto triangleNormal(int t, bool normalize=true) const
SYS_FORCE_INLINE int dstVertex(int h) const
SYS_FORCE_INLINE auto dualHedgeLength(int h) const
SYS_FORCE_INLINE auto triangleBarycenter(int t) const
SYS_FORCE_INLINE auto pointDualArea(int p) const
UT_Vector3T< T > rotate(const UT_Vector3T< T > &) const
#define SYSmin(a, b)
Definition: SYS_Math.h:1448
SYS_FORCE_INLINE auto hedgeVector(int h) const
int numInteriorHedges() const
SYS_FORCE_INLINE int hedgeBoundary(int h) const
Declare prior to use.
SYS_FORCE_INLINE int trianglePoint(int t, int j=0) const
SYS_FORCE_INLINE auto hedgeCotan(int h) const
void UTparallelFor(const Range &range, const Body &body, const int subscribe_ratio=2, const int min_grain_size=1)
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
GLdouble GLdouble t
Definition: glew.h:1398
GLuint res
Definition: glew.h:11507
SYS_FORCE_INLINE int srcPoint(int h) const
SYS_FORCE_INLINE GA_Offset dstPoint(T &iface, GEO_Hedge h)
Definition: GEO_Hedge.h:244
SYS_FORCE_INLINE bool isBoundaryHedge(T &iface, GEO_Hedge h)
Definition: GEO_Hedge.h:309