1+ #include " geometrycentral/surface/marching_triangles.h"
2+
3+ namespace geometrycentral {
4+ namespace surface {
5+
6+ namespace {
7+
8+ std::vector<std::vector<std::array<size_t , 2 >>>
9+ getCurveComponents (SurfaceMesh& mesh, const std::vector<SurfacePoint>& curveNodes,
10+ const std::vector<std::array<size_t , 2 >>& curveEdges) {
11+
12+ // Note: This function assumes that `curveNodes` has been de-duplicated.
13+ std::vector<std::array<size_t , 2 >> edgesToAdd = curveEdges;
14+ std::vector<std::vector<std::array<size_t , 2 >>> curves;
15+ size_t nSegs = curveEdges.size ();
16+ while (edgesToAdd.size () > 0 ) {
17+ std::array<size_t , 2 > startSeg = edgesToAdd.back ();
18+ edgesToAdd.pop_back ();
19+ curves.emplace_back ();
20+ std::vector<std::array<size_t , 2 >>& currCurve = curves.back ();
21+ currCurve.push_back (startSeg);
22+
23+ // Add segs to the front end until we can't.
24+ std::array<size_t , 2 > currSeg = startSeg;
25+ while (true ) {
26+ const SurfacePoint& front = curveNodes[currSeg[1 ]];
27+ bool didWeFindOne = false ;
28+ for (size_t i = 0 ; i < edgesToAdd.size (); i++) {
29+ std::array<size_t , 2 > otherSeg = edgesToAdd[i];
30+ if (curveNodes[otherSeg[0 ]] == front) {
31+ currSeg = otherSeg;
32+ currCurve.push_back (otherSeg);
33+ edgesToAdd.erase (edgesToAdd.begin () + i);
34+ didWeFindOne = true ;
35+ break ;
36+ }
37+ }
38+ if (!didWeFindOne) break ;
39+ }
40+
41+ // Add segs to the back end until we can't.
42+ currSeg = startSeg;
43+ while (true ) {
44+ const SurfacePoint& back = curveNodes[currSeg[0 ]];
45+ bool didWeFindOne = false ;
46+ for (size_t i = 0 ; i < edgesToAdd.size (); i++) {
47+ std::array<size_t , 2 > otherSeg = edgesToAdd[i];
48+ if (curveNodes[otherSeg[1 ]] == back) {
49+ currSeg = otherSeg;
50+ currCurve.insert (currCurve.begin (), otherSeg);
51+ edgesToAdd.erase (edgesToAdd.begin () + i);
52+ didWeFindOne = true ;
53+ break ;
54+ }
55+ }
56+ if (!didWeFindOne) break ;
57+ }
58+ }
59+ return curves;
60+ }
61+ } // namespace
62+
63+ std::vector<std::vector<SurfacePoint>> marchingTriangles (IntrinsicGeometryInterface& geom, const VertexData<double >& u,
64+ double isoval) {
65+
66+ SurfaceMesh& mesh = *u.getMesh ();
67+ std::vector<SurfacePoint> nodes;
68+ std::vector<std::array<size_t , 2 >> edges;
69+ for (Face f : mesh.faces ()) {
70+ std::vector<SurfacePoint> hits;
71+ BarycentricVector gradient (f);
72+ for (Halfedge he : f.adjacentHalfedges ()) {
73+ // Record edge crossings
74+ Edge e = he.edge ();
75+ Vertex v0 = e.firstVertex ();
76+ Vertex v1 = e.secondVertex ();
77+ double u0 = u[v0];
78+ double u1 = u[v1];
79+ double lB = std::min (u0, u1);
80+ double uB = std::max (u0, u1);
81+ if (lB == uB && lB == isoval) {
82+ hits.clear ();
83+ hits = {SurfacePoint (v0), SurfacePoint (v1)};
84+ break ;
85+ }
86+ double t = (isoval - lB) / (uB - lB);
87+ if (u0 > u1) t = 1 . - t;
88+ if (t <= 1 . && t >= 0 .) hits.emplace_back (e, t);
89+ // Compute gradient of the scalar function.
90+ BarycentricVector heVec (he.next (), f);
91+ BarycentricVector ePerp = heVec.rotate90 (geom);
92+ gradient += ePerp * u[he.vertex ()];
93+ }
94+ if (hits.size () != 2 ) continue ;
95+
96+ // Orient segments so that smaller values are always on the "inside" of the curve.
97+ std::array<size_t , 2 > seg;
98+ for (int i = 0 ; i < 2 ; i++) {
99+ auto iter = std::find (nodes.begin (), nodes.end (), hits[i]);
100+ if (iter != nodes.end ()) {
101+ seg[i] = iter - nodes.begin ();
102+ } else {
103+ nodes.push_back (hits[i]);
104+ seg[i] = nodes.size () - 1 ;
105+ }
106+ }
107+ BarycentricVector segTangent (hits[0 ], hits[1 ]);
108+ BarycentricVector segNormal = -segTangent.inFace (f).rotate90 (geom);
109+ if (dot (geom, segNormal, gradient) > 0 .) {
110+ edges.push_back (seg);
111+ } else {
112+ edges.push_back ({seg[1 ], seg[0 ]});
113+ }
114+ }
115+ std::vector<std::vector<std::array<size_t , 2 >>> components = getCurveComponents (mesh, nodes, edges);
116+ size_t nCurves = components.size ();
117+ std::vector<std::vector<SurfacePoint>> curves (nCurves);
118+ for (size_t i = 0 ; i < nCurves; i++) {
119+ size_t nEdges = components[i].size ();
120+ for (size_t j = 0 ; j < nEdges; j++) {
121+ curves[i].push_back (nodes[components[i][j][0 ]]);
122+ }
123+ curves[i].push_back (nodes[components[i][nEdges - 1 ][1 ]]);
124+ }
125+ return curves;
126+ }
127+
128+ } // namespace surface
129+ } // namespace geometrycentral
0 commit comments