-
-
Save simonlindholm/050cfd03f7def7e8714123953a16432e to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/** | |
* Author: Simon Lindholm | |
* Date: 2020-10-17 | |
* Source: http://www.inf.fu-berlin.de/lehre/SS13/AG/randinc3d.pdf | |
* Description: Computes all faces of the 3-dimension hull of a point set. | |
* \begin{enumerate} | |
* \item \texttt{verts} gives vertices of each triangle, in ccw order seen from outside. | |
* \item \texttt{edges[i]} is the adjacent triangle when crossing the (vertex i+1, vertex i+2) edge. | |
* \end{enumerate} | |
* Time: expected O(N log N) | |
* Status: TODO stress-test and SPOJ CH3D | |
* | |
* Details: | |
* The basic algorithm is the following: | |
* - Shuffle the input points randomly, making sure the first four points are | |
* in general position. Set up a tetrahedron from them. | |
* - For all other points, try to find an arbitrary triangle that it lies | |
* outside of ("conflicts with"), i.e. which should be removed when the point | |
* is added. Keep these in list 'dep', plus a reverse list 'rdep' from | |
* triangle to points. | |
* - Loop over the other points, adding them one by one if they have a conflict | |
* triangle (points without conflict triangles are not in the hull). A point | |
* can have multiple conflict triangles, and we can find them all by | |
* performing a BFS from dep[i]. Create new triangles from the new point and | |
* the border of the BFS'd region. (At this point, we have edges pointing | |
* from removed triangle -> new triangle -> border triangle -> new triangle.) | |
* This part takes O(n) time in total, in expectation. | |
* - For all removed triangles, recompute conflict triangles for the points in | |
* rdep[]. This is done by a BFS over the set of removed triangles plus the | |
* newly added ones, looking at only conflict triangles, starting from the | |
* known conflict triangle. This part takes O(n log n) time. | |
* - Finally, reindex back into input coordinates. | |
* | |
* We keep the following state: | |
* - 'pi': a random permutation of 0..N, such that ps[pi[0..4]] (used as the | |
* starting tetrahedron) are in general position. | |
* - 'pp': pp[i] = ps[pi[i]]. This is an optimization to reduce indirection. | |
* - 'tris': a buffer of triangle hull faces, some of which are live. Unused | |
* triangles are kept in a freelist `free'. In each step we may need up to 2N | |
* triangles: V-E+F=2, E=3/2*F, F = 2V-4 < 2N. Since both removed and new | |
* triangles are live at once, we may need twice that temporarily (or perhaps | |
* only 1.5x? but 4N is cheap enough). Hence, we give it size T = 4*N. The | |
* triangle with index 0 is never used, to get nice sentinel for dep[], and | |
* to help the partial_sum-based reindexing at the end. | |
* - 'free': a stack of unused triangles, with logical size 'fend'. It does | |
* contain 0, to make the loop at the end work, but it's at the bottom of the | |
* stack and will never be used. | |
* - 'dep', 'rdep': explained above. | |
* - 'X': when removing the triangles bounded by the cycle Z, and introducing | |
* new ones for all edges on the cycle, X[vert] = tri1 ^ tri2 where vert is a | |
* vertex of the cycle, and tri1 and tri2 are the two newly added triangles. | |
* This helps the triangles find each other. At the start of each loop | |
* iteration, X will be all zero. | |
* - 'seen': when performing BFS's/DFS's, this keeps track of which triangles | |
* have been seen; ones unseen have seen[i] < time. time increments for each | |
* search. | |
* - 'q': queue for the initial BFS. Each triangle is added just once into the | |
* queue, at the same time as marking it seen. 'qe' marks the end of the | |
* queue. The queue is later traversed to recompute conflict triangles, which | |
* is why we cannot reuse the same queue for the DFS's, and why we cannot use | |
* a DFS for the first traversal. | |
* - 's': stack for the subsequent DFS's. 'e' marks the end. | |
* - 'ren': mapping from sparse triangle indices to dense ones that we return. | |
*/ | |
#pragma once | |
struct Tri { array<int, 3> verts, edges; }; | |
vector<Tri> hull3d(const vector<Point3D<double>>& ps) { | |
int N = sz(ps), T = 4 * N, time = 0, fend = 1, p = 2, tri = 4; | |
auto pp = ps; | |
vi pi(N), X(N), dep(N), seen(T), free(T), q(T), s(T), ren(T, 1); | |
vector<Tri> tris(T), v; | |
rep(i,0,4) rep(j,0,3) | |
++(tris[i+1].edges[j] = tris[i+1].verts[j] = i ^ (j + 1)); | |
iota(all(pi), 0); | |
while (!(ps[p] - ps[0]).cross(ps[1] - ps[0]).dist2()) p++; // TODO: don't use dist2() | |
pp[2] = pp[p]; swap(pi[p++], pi[2]); | |
#define PT(i) (pp[tris[tri].verts[i]] - pp[p]) | |
auto outside = [&](int tri, int p) { | |
return PT(0).cross(PT(1)).dot(PT(2)) < 0; // TODO: int128? | |
}; | |
while (!PT(0).cross(PT(1)).dot(PT(2))) p++, assert(p < N); // TODO: int128? | |
swap(pi[3], pi[p]); | |
if (outside(tri, p)) swap(pi[0], pi[1]); | |
random_shuffle(4 + all(pi)); | |
rep(i,0,N) pp[i] = ps[pi[i]]; | |
rep(i,5,T) free[fend++] = i; | |
vector<vi> rdep(T); | |
rep(i,0,4) rep(j,4,N) if (!dep[j] && outside(i, j)) | |
dep[j] = i, rdep[i].push_back(j); | |
rep(j,4,N) if (dep[j]) { | |
for (int x : X) assert(!x); | |
int t = q[0] = dep[j], qe = 1, origInd = fend; | |
seen[t] = ++time; | |
assert(outside(q[0], j)); | |
rep(qi,0,qe) rep(ei,0,3) { | |
int t = q[qi], &t2 = tris[t].edges[ei]; | |
if (seen[t2] == time) continue; | |
if (outside(t2, j)) { seen[q[qe++] = t2] = time; continue; } | |
int newTri = free[--fend]; | |
rdep[newTri].clear(); | |
Tri& tr = tris[newTri]; | |
tr.verts[0] = j; | |
tr.edges[0] = t2; | |
rep(i,1,3) tr.verts[i] = tris[t].verts[(ei+i) % 3]; | |
rep(i,1,3) X[tr.verts[i]] ^= newTri; | |
replace(all(tris[t2].edges), t, newTri); | |
t2 = newTri; | |
} | |
rep(i,fend,origInd) { | |
int t = free[i], &x = X[tris[t].verts[1]]; | |
tris[tris[t].edges[2] = t ^ x].edges[1] = t; | |
x = 0; | |
} | |
rep(qi,0,qe) free[fend++] = q[qi]; | |
rep(qi,0,qe) for (int x : rdep[q[qi]]) if (x != j) { | |
int e = 1; s[0] = dep[x]; | |
assert(outside(dep[x], x)); | |
seen[dep[x]] = ++time; | |
dep[x] = 0; | |
while (e-- && tris[s[e]].verts[0] != j) // TODO: could also check "isOld[s[e]]" | |
for (int t : tris[s[e]].edges) | |
if (seen[t] < time && outside(t, x)) | |
s[e++] = t, seen[t] = time; | |
if (e != -1) rdep[dep[x] = s[e]].push_back(x); | |
} | |
} | |
rep(i,0,fend) ren[free[i]] = 0; | |
partial_sum(all(ren), ren.begin()); | |
v.resize(ren.back()); | |
rep(i,1,T) if (ren[i] != ren[i - 1]) { | |
for (int& x : tris[i].edges) x = ren[x - 1]; | |
for (int& x : tris[i].verts) x = pi[x]; | |
v[ren[i - 1]] = tris[i]; | |
} | |
return v; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment