diff options
Diffstat (limited to 'geometry')
| -rw-r--r-- | geometry/delaunay.cpp | 124 | ||||
| -rw-r--r-- | geometry/geometry.tex | 34 | ||||
| -rw-r--r-- | geometry/hpi.cpp | 68 | ||||
| -rw-r--r-- | geometry/linesAndSegments.cpp | 3 | ||||
| -rw-r--r-- | geometry/polygon.cpp | 6 | ||||
| -rw-r--r-- | geometry/triangle.cpp | 15 |
6 files changed, 229 insertions, 21 deletions
diff --git a/geometry/delaunay.cpp b/geometry/delaunay.cpp new file mode 100644 index 0000000..1008b39 --- /dev/null +++ b/geometry/delaunay.cpp @@ -0,0 +1,124 @@ +using lll = __int128; +using pt = complex<ll>; + +constexpr pt INF_PT = pt(1e18, 1e18); + +bool circ(pt p, pt a, pt b, pt c) {// p in circle(A,B,C) + return imag((c-b)*conj(p-c)*(a-p)*conj(b-a)) < 0; +} + +struct QuadEdge { + QuadEdge* rot = nullptr; + QuadEdge* onext = nullptr; + pt orig = INF_PT; + bool used = false; + QuadEdge* rev() const {return rot->rot;} + QuadEdge* lnext() const {return rot->rev()->onext->rot;} + QuadEdge* oprev() const {return rot->onext->rot;} + pt dest() const {return rev()->orig;} +}; + +deque<QuadEdge> edgeData; + +QuadEdge* makeEdge(pt from, pt to) { + for (int j : {0,1,2,3}) edgeData.push_back({}); + auto e = edgeData.end() - 4; + for (int j : {0,1,2,3}) e[j].onext = e[j^3].rot = &e[j^(j>>1)]; + e[0].orig = from; + e[1].orig = to; + return &e[0]; +} + +void splice(QuadEdge* a, QuadEdge* b) { + swap(a->onext->rot->onext, b->onext->rot->onext); + swap(a->onext, b->onext); +} + +QuadEdge* connect(QuadEdge* a, QuadEdge* b) { + QuadEdge* e = makeEdge(a->dest(), b->orig); + splice(e, a->lnext()); + splice(e->rev(), b); + return e; +} + +bool valid(QuadEdge* e, QuadEdge* base) { + return cross(e->dest(), base->orig, base->dest()) < 0; +} + +template<bool ccw> +QuadEdge* deleteAll(QuadEdge* e, QuadEdge* base) { + if (valid(e, base)) { + while (circ(base->dest(), base->orig, e->dest(), (ccw ? e->onext : e->oprev())->dest())) { + QuadEdge* t = ccw ? e->onext : e->oprev(); + splice(e, e->oprev()); + splice(e->rev(), e->rev()->oprev()); + e = t; + }} + return e; +} + +template<typename IT> +pair<QuadEdge*, QuadEdge*> rec(IT l, IT r) { + int n = distance(l, r); + if (n <= 3) { + QuadEdge* a = makeEdge(l[0], l[1]); + if (n == 2) return {a, a->rev()}; + QuadEdge* b = makeEdge(l[1], l[2]); + splice(a->rev(), b); + int side = cross(l[0], l[1], l[2]); + QuadEdge* c = nullptr; + if (side != 0) c = connect(b, a); + if (side >= 0) return {a, b->rev()}; + else return {c->rev(), c}; + } + auto m = l + (n / 2); + auto [ldo, ldi] = rec(l, m); + auto [rdi, rdo] = rec(m, r); + while (true) { + if (cross(rdi->orig, ldi->orig, ldi->dest()) > 0) { + ldi = ldi->lnext(); + } else if (cross(ldi->orig, rdi->orig, rdi->dest()) < 0) { + rdi = rdi->rev()->onext; + } else break; + } + QuadEdge* base = connect(rdi->rev(), ldi); + if (ldi->orig == ldo->orig) ldo = base->rev(); + if (rdi->orig == rdo->orig) rdo = base; + while (true) { + QuadEdge* lcand = deleteAll<true>(base->rev()->onext, base); + QuadEdge* rcand = deleteAll<false>(base->oprev(), base); + if (!valid(lcand, base) && !valid(rcand, base)) break; + if (!valid(lcand, base) || (valid(rcand, base) && + circ(lcand->dest(), lcand->orig, rcand->orig, rcand->dest()))) { + base = connect(rcand, base->rev()); + } else { + base = connect(base->rev(), lcand->rev()); + }} + return {ldo, rdo}; +} + +vector<pt> delaunay(vector<pt> pts) { + if (sz(pts) <= 2) return {}; + sort(all(pts), [](const pt& a, const pt& b) { + if (real(a) != real(b)) return real(a) < real(b); + return imag(a) < imag(b); + }); + QuadEdge* r = rec(all(pts)).first; + vector<QuadEdge*> edges = {r}; + while (cross(r->onext->dest(), r->dest(), r->orig) < 0) r = r->onext; + auto add = [&](QuadEdge* e){ + QuadEdge* cur = e; + do { + cur->used = true; + pts.push_back(cur->orig); + edges.push_back(cur->rev()); + cur = cur->lnext(); + } while (cur != e); + }; + add(r); + pts.clear(); + for (int i = 0; i < sz(edges); i++) { + if (!edges[i]->used) add(edges[i]); + } + return pts; +} diff --git a/geometry/geometry.tex b/geometry/geometry.tex index a027de4..d92bad1 100644 --- a/geometry/geometry.tex +++ b/geometry/geometry.tex @@ -7,14 +7,6 @@ \sourcecode{geometry/closestPair.cpp} \end{algorithm} -\begin{algorithm}{Rotating calipers} - \begin{methods} - \method{antipodalPoints}{berechnet antipodale Punkte}{n} - \end{methods} - \textbf{WICHTIG:} Punkte müssen gegen den Uhrzeigersinn Sortiert sein und konvexes Polygon bilden! - \sourcecode{geometry/antipodalPoints.cpp} -\end{algorithm} - \begin{algorithm}{Konvexe Hülle} \begin{methods} \method{convexHull}{berechnet Konvexehülle}{n\*\log(n)} @@ -27,23 +19,43 @@ \sourcecode{geometry/convexHull.cpp} \end{algorithm} +\begin{algorithm}{Rotating calipers} + \begin{methods} + \method{antipodalPoints}{berechnet antipodale Punkte}{n} + \end{methods} + \textbf{WICHTIG:} Punkte müssen gegen den Uhrzeigersinn Sortiert sein und konvexes Polygon bilden! + \sourcecode{geometry/antipodalPoints.cpp} +\end{algorithm} + \subsection{Formeln~~--~\texttt{std::complex}} \sourcecode{geometry/formulars.cpp} \sourcecode{geometry/linesAndSegments.cpp} +\sourcecode{geometry/sortAround.cpp} \input{geometry/triangle} \sourcecode{geometry/triangle.cpp} \sourcecode{geometry/polygon.cpp} -\sourcecode{geometry/sortAround.cpp} \sourcecode{geometry/circle.cpp} \subsection{Formeln - 3D} \sourcecode{geometry/formulars3d.cpp} \optional{ -\subsection{3D-Kugeln} -\sourcecode{geometry/spheres.cpp} + \subsection{3D-Kugeln} + \sourcecode{geometry/spheres.cpp} } +\begin{algorithm}{Half-plane intersection} + \sourcecode{geometry/hpi.cpp} +\end{algorithm} + +\begin{algorithm}[optional]{Delaunay Triangulierung} + \begin{methods} + \method{delaunay}{berechnet Triangulierung}{n\*\log(n)} + \end{methods} + \textbf{WICHTIG:} Wenn alle Punkte kollinear sind gibt es keine Traingulierung! Wenn 4 Punkte auf einem Kreis liegen ist die Triangulierung nicht eindeutig. + \sourcecode{geometry/delaunay.cpp} +\end{algorithm} + \optional{ \subsection{Geraden} \sourcecode{geometry/lines.cpp} diff --git a/geometry/hpi.cpp b/geometry/hpi.cpp new file mode 100644 index 0000000..3509e0e --- /dev/null +++ b/geometry/hpi.cpp @@ -0,0 +1,68 @@ +constexpr ll inf = 0x1FFF'FFFF'FFFF'FFFF;//THIS CODE IS WIP + +bool left(pt p) {return real(p) < 0 || + (real(p) == 0 && imag(p) < 0);} +struct hp { + pt from, to; + + hp(pt a, pt b) : from(a), to(b) {} + hp(pt dummy) : hp(dummy, dummy) {} + + bool dummy() const {return from == to;} + pt dir() const {return dummy() ? to : to - from;} + bool operator<(const hp& o) const { + if (left(dir()) != left(o.dir())) + return left(dir()) > left(o.dir()); + return cross(dir(), o.dir()) > 0; + } + + using lll = __int128; + using ptl = complex<lll>; + ptl mul(lll m, ptl p) const {return m*p;}//ensure 128bit + + bool check(const hp& a, const hp& b) const { + if (dummy() || b.dummy()) return false; + if (a.dummy()) { + ll ort = sgn(cross(b.dir(), dir())); + if (ort == 0) return cross(from, to, a.from) < 0; + return cross(b.dir(), a.dir()) * ort > 0; + } + ll y = cross(a.dir(), b.dir()); + ll z = cross(b.from - a.from, b.dir()); + ptl i = mul(y, a.from) + mul(z, a.dir()); //intersect a and b + // check if i is outside/right of x + return imag(conj(mul(sgn(y),dir()))*(i-mul(y,from))) < 0; + } +}; + +constexpr ll lim = 2e9+7; + +deque<hp> intersect(vector<hp> hps) { + hps.push_back(hp(pt{lim+1,-1})); + hps.push_back(hp(pt{lim+1,1})); + sort(all(hps)); + + deque<hp> dq = {hp(pt{-lim, 1})}; + for (auto x : hps) { + while (sz(dq) > 1 && x.check(dq.end()[-1], dq.end()[-2])) + dq.pop_back(); + while (sz(dq) > 1 && x.check(dq[0], dq[1])) + dq.pop_front(); + + if (cross(x.dir(), dq.back().dir()) == 0) { + if (dot(x.dir(), dq.back().dir()) < 0) return {}; + if (cross(x.from, x.to, dq.back().from) < 0) + dq.pop_back(); + else continue; + } + dq.push_back(x); + } + + while (sz(dq) > 2 && dq[0].check(dq.end()[-1], dq.end()[-2])) + dq.pop_back(); + while (sz(dq) > 2 && dq.end()[-1].check(dq[0], dq[1])) + dq.pop_front(); + + if (sz(dq) < 3) return {}; + return dq; +} diff --git a/geometry/linesAndSegments.cpp b/geometry/linesAndSegments.cpp index ba6c468..c86b331 100644 --- a/geometry/linesAndSegments.cpp +++ b/geometry/linesAndSegments.cpp @@ -17,8 +17,7 @@ vector<pt> lineSegmentIntersection(pt p0, pt p1, pt p2, pt p3) { double b = cross(p2 - p0, p3 - p2); double c = cross(p1 - p0, p0 - p2); if (a < 0) {a = -a; b = -b; c = -c;} - if (b < -EPS || a-b < -EPS || - c < -EPS || a-c < -EPS) return {}; + if (b < -EPS || b-a > EPS || c < -EPS || c-a > EPS) return {}; if (a > EPS) return {p0 + b/a*(p1 - p0)}; vector<pt> result; auto insertUnique = [&](pt p) { diff --git a/geometry/polygon.cpp b/geometry/polygon.cpp index f562b7a..409f765 100644 --- a/geometry/polygon.cpp +++ b/geometry/polygon.cpp @@ -10,13 +10,13 @@ double area(const vector<pt>& poly) { //poly[0] == poly.back() // Anzahl drehungen einer Polyline um einen Punkt // p nicht auf rand und poly[0] == poly.back() // res != 0 or (res & 1) != 0 um inside zu prüfen bei -// selbstschneidenden polygonen (definitions sache) +// selbstschneidenden Polygonen (definitions Sache) ll windingNumber(pt p, const vector<pt>& poly) { ll res = 0; for (int i = 0; i + 1 < sz(poly); i++) { pt a = poly[i], b = poly[i + 1]; - if (real(a) > real(b)) swap(a,b); - if (real(a) <= real(p) &&real(p) < real(b) && + if (real(a) > real(b)) swap(a, b); + if (real(a) <= real(p) && real(p) < real(b) && cross(p, a, b) < 0) { res += orientation(p, poly[i], poly[i + 1]); }} diff --git a/geometry/triangle.cpp b/geometry/triangle.cpp index a00eb56..ebbbc88 100644 --- a/geometry/triangle.cpp +++ b/geometry/triangle.cpp @@ -28,12 +28,17 @@ pt outCenter(pt a, pt b, pt c) { c*conj(c)*conj(a-b)) / d; } +// 1 => p außerhalb Kreis durch a,b,c +// 0 => p auf Kreis durch a,b,c +// -1 => p im Kreis durch a,b,c +int insideOutCenter(pt a, pt b, pt c, pt p) {// braucht lll + return sgn(imag((c-b)*conj(p-c)*(a-p)*conj(b-a))); +} + // Sind die Dreiecke a1, b1, c1, and a2, b2, c2 ähnlich? // Erste Zeile testet Ähnlichkeit mit gleicher Orientierung, // zweite Zeile testet Ähnlichkeit mit verschiedener Orientierung -bool similar (pt a1, pt b1, pt c1, pt a2, pt b2, pt c2) { - return ((b2-a2) * (c1-a1) == (b1-a1) * (c2-a2) || - (b2-a2) * (conj(c1)-conj(a1)) == (conj(b1)-conj(a1)) - * (c2-a2) - ); +bool similar(pt a1, pt b1, pt c1, pt a2, pt b2, pt c2) { + return (b2-a2) * (c1-a1) == (b1-a1) * (c2-a2) || + (b2-a2) * conj(c1-a1) == conj(b1-a1) * (c2-a2); } |