diff options
Diffstat (limited to 'geometry')
| -rw-r--r-- | geometry/antipodalPoints.cpp | 13 | ||||
| -rw-r--r-- | geometry/circle.cpp | 33 | ||||
| -rw-r--r-- | geometry/closestPair.cpp | 39 | ||||
| -rw-r--r-- | geometry/convexHull.cpp | 39 | ||||
| -rw-r--r-- | geometry/formulars.cpp | 174 | ||||
| -rw-r--r-- | geometry/formulars3d.cpp | 53 | ||||
| -rw-r--r-- | geometry/geometry.tex | 50 | ||||
| -rw-r--r-- | geometry/lines.cpp | 29 | ||||
| -rw-r--r-- | geometry/linesAndSegments.cpp | 90 | ||||
| -rw-r--r-- | geometry/polygon.cpp | 33 | ||||
| -rw-r--r-- | geometry/segmentIntersection.cpp | 63 | ||||
| -rw-r--r-- | geometry/sortAround.cpp | 10 | ||||
| -rw-r--r-- | geometry/spheres.cpp | 29 | ||||
| -rw-r--r-- | geometry/triangle.cpp | 40 |
14 files changed, 483 insertions, 212 deletions
diff --git a/geometry/antipodalPoints.cpp b/geometry/antipodalPoints.cpp new file mode 100644 index 0000000..c1921cc --- /dev/null +++ b/geometry/antipodalPoints.cpp @@ -0,0 +1,13 @@ +vector<pair<int, int>> antipodalPoints(vector<pt>& h) { + vector<pair<int, int>> result; + int n = (int)h.size(); + if (n < 2) return result; + for (int i = 0, j = 1; i < j; i++) { + while (true) { + result.push_back({i, j}); + if (cross(h[(i + 1) % n] - h[i], + h[(j + 1) % n] - h[j]) <= 0) break; + j = (j + 1) % n; + }} + return result; +} diff --git a/geometry/circle.cpp b/geometry/circle.cpp new file mode 100644 index 0000000..f49ea27 --- /dev/null +++ b/geometry/circle.cpp @@ -0,0 +1,33 @@ +// berechnet die Schnittpunkte von zwei kreisen +// (Kreise dürfen nicht gleich sein!) +vector<pt> circleIntersection(pt c1, double r1, + pt c2, double r2) { + double d = abs(c1 - c2); + if (d < abs(r1 - r2) || d > abs(r1 + r2)) return {}; + double a = (r1 * r1 - r2 * r2 + d * d) / (2 * d); + pt p = (c2 - c1) * a / d + c1; + if (d == abs(r1 - r2) || d == abs(r1 + r2)) return {p}; + double h = sqrt(r1 * r1 - a * a); + return {p + pt{0, 1} * (c2 - c1) * h / d, + p - pt{0, 1} * (c2 - c1) * h / d}; +} + +// berechnet die Schnittpunkte zwischen +// einem Kreis(Kugel) und einer Grade 2d und 3d +vector<pt> circleRayIntersection(pt center, double r, + pt orig, pt dir) { + vector<pt> result; + double a = dot(dir, dir); + double b = 2 * dot(dir, orig - center); + double c = dot(orig - center, orig - center) - r * r; + double discr = b * b - 4 * a * c; + if (discr >= 0) { + //t in [0, 1] => schnitt mit segment [orig, orig + dir] + double t1 = -(b + sqrt(discr)) / (2 * a); + double t2 = -(b - sqrt(discr)) / (2 * a); + if (t1 >= 0) result.push_back(t1 * dir + orig); + if (t2 >= 0 && abs(t1 - t2) > EPS) { + result.push_back(t2 * dir + orig); + }} + return result; +} diff --git a/geometry/closestPair.cpp b/geometry/closestPair.cpp index cac544b..a7662b6 100644 --- a/geometry/closestPair.cpp +++ b/geometry/closestPair.cpp @@ -1,35 +1,40 @@ double squaredDist(pt a, pt b) { - return (a.fst-b.fst) * (a.fst-b.fst) + (a.snd-b.snd) * (a.snd-b.snd); + return real(conj(a-b) * (a-b)); } bool compY(pt a, pt b) { - if (a.snd == b.snd) return a.fst < b.fst; - return a.snd < b.snd; + return (imag(a) == imag(b)) ? real(a) < real(b) + : imag(a) < imag(b); } -// points.size() > 1 und alle Punkte müssen verschieden sein! -double shortestDist(vector<pt> &points) { +bool compX(pt a, pt b) { + return (real(a) == real(b)) ? imag(a) < imag(b) + : real(a) < real(b); +} + +double shortestDist(vector<pt>& points) { // points.size() > 1 set<pt, bool(*)(pt, pt)> status(compY); - sort(points.begin(), points.end()); - double opt = 1e30, sqrtOpt = 1e15; + sort(points.begin(), points.end(), compX); + double opt = 1.0/0.0, sqrtOpt = 1.0/0.0; auto left = points.begin(), right = points.begin(); status.insert(*right); right++; - + while (right != points.end()) { - if (fabs(left->fst - right->fst) >= sqrtOpt) { - status.erase(*(left++)); + if (left != right && + abs(real(*left - *right)) >= sqrtOpt) { + status.erase(*left); + left++; } else { - auto lower = status.lower_bound(pt(-1e20, right->snd - sqrtOpt)); - auto upper = status.upper_bound(pt(-1e20, right->snd + sqrtOpt)); - while (lower != upper) { + auto lower = status.lower_bound({-1.0/0.0, imag(*right) - sqrtOpt}); + auto upper = status.upper_bound({-1.0/0.0, imag(*right) + sqrtOpt}); + for (;lower != upper; lower++) { double cand = squaredDist(*right, *lower); if (cand < opt) { opt = cand; sqrtOpt = sqrt(opt); - } - ++lower; - } - status.insert(*(right++)); + }} + status.insert(*right); + right++; }} return sqrtOpt; } diff --git a/geometry/convexHull.cpp b/geometry/convexHull.cpp index 9c14c45..e988977 100644 --- a/geometry/convexHull.cpp +++ b/geometry/convexHull.cpp @@ -1,24 +1,19 @@ -// Laufzeit: O(n*log(n)) -ll cross(const pt p, const pt a, const pt b) { - return (a.x - p.x) * (b.y - p.y) - (a.y - p.y) * (b.x - p.x); -} - -// Punkte auf der konvexen Hülle, gegen den Uhrzeigersinn sortiert. -// Kollineare Punkte nicht enthalten, entferne dafür "=" im CCW-Test. -// Achtung: Der erste und letzte Punkt im Ergebnis sind gleich. -// Achtung: Alle Punkte müssen verschieden sein. vector<pt> convexHull(vector<pt> p){ - int n = p.size(), k = 0; - vector<pt> h(2 * n); - sort(p.begin(), p.end()); - for (int i = 0; i < n; i++) { // Untere Hülle. - while (k >= 2 && cross(h[k - 2], h[k - 1], p[i]) <= 0.0) k--; - h[k++] = p[i]; - } - for (int i = n - 2, t = k + 1; i >= 0; i--) { // Obere Hülle. - while (k >= t && cross(h[k - 2], h[k - 1], p[i]) <= 0.0) k--; - h[k++] = p[i]; - } - h.resize(k); - return h; + sort(p.begin(), p.end(), [](const pt& a, const pt& b){ + return real(a) == real(b) ? imag(a) < imag(b) + : real(a) < real(b); + }); + p.erase(unique(p.begin(), p.end()), p.end()); + int k = 0; + vector<pt> h(2 * sz(p)); + for (int i = 0; i < sz(p); i++) {// Untere Hülle. + while (k > 1 && cross(h[k-2], h[k-1], p[i]) <= 0) k--; + h[k++] = p[i]; + } + for (int i = sz(p)-2, t = k; i >= 0; i--) {// Obere Hülle. + while (k > t && cross(h[k-2], h[k-1], p[i]) <= 0) k--; + h[k++] = p[i]; + } + h.resize(k); + return h; } diff --git a/geometry/formulars.cpp b/geometry/formulars.cpp index ed44b8b..43affbb 100644 --- a/geometry/formulars.cpp +++ b/geometry/formulars.cpp @@ -1,165 +1,37 @@ -// Komplexe Zahlen als Darstellung für Punkte. Wenn immer möglich -// complex<int> verwenden. Funktionen wie abs() geben dann int zurück. -typedef complex<double> pt; +// Komplexe Zahlen als Darstellung für Punkte. +// Wenn immer möglich complex<int> verwenden. +// Funktionen wie abs() geben dann int zurück. +using pt = complex<double>; -// Winkel zwischen Punkt und x-Achse in [0, 2 * PI). -double angle = arg(a); +// PIL < PI < PIU +constexpr double PIU = acos(-1.0l); +constexpr double PIL = PIU-2e-19l; -// Punkt rotiert um Winkel theta. -pt a_rotated = a * exp(pt(0, theta)); +// Winkel zwischen Punkt und x-Achse in [-PI, PI]. +double angle(pt a) {return arg(a);} -// Mittelpunkt des Dreiecks abc. -pt centroid = (a + b + c) / 3.0; +// rotiert Punkt im Uhrzeigersinn um den Ursprung. +pt rotate(pt a, double theta) {return a * exp(pt(0.0, theta));} // Skalarprodukt. -double dot(pt a, pt b) { return real(conj(a) * b); } +double dot(pt a, pt b) {return real(conj(a) * b);} -// Kreuzprodukt, 0, falls kollinear. -double cross(pt a, pt b) { return imag(conj(a) * b); } - -// Flächeninhalt eines Dreicks bei bekannten Eckpunkten. -double areaOfTriangle(pt a, pt b, pt c) { - return abs(cross(b - a, c - a)) / 2.0; -} +// abs()^2.(pre c++20) +double norm(pt a) {return dot(a, a);} -// Flächeninhalt eines Dreiecks bei bekannten Seitenlängen. -double areaOfTriangle(double a, double b, double c) { - double s = (a + b + c) / 2; - return sqrt(s * (s-a) * (s-b) * (s-c)); -} - -// Sind die Dreiecke a1, b1, c1, and a2, b2, c2 ähnlich? -// Erste Zeile testet Ähnlichkeit mit gleicher Orientierung, -// zweite Zeile testet Ähnlichkeit mit unterschiedlicher 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) - ); -} +// Kreuzprodukt, 0, falls kollinear. +double cross(pt a, pt b) {return imag(conj(a) * b);} +double cross(pt p, pt a, pt b) {return cross(a - p, b - p);} -// -1 => gegen den Uhrzeigersinn, 0 => kolliniear, 1 => im Uhrzeigersinn. -// Einschränken der Rückgabe auf [-1,1] ist sicherer gegen Overflows. -double orientation(pt a, pt b, pt c) { +// -1 => gegen den Uhrzeigersinn +// 0 => kolliniear +// 1 => im Uhrzeigersinn. +int orientation(pt a, pt b, pt c) { double orien = cross(b - a, c - a); - if (abs(orien) < EPSILON) return 0; // Braucht großes EPSILON: ~1e-6 - return orien < 0 ? -1 : 1; -} - -// Test auf Streckenschnitt zwischen a-b und c-d. -bool lineSegmentIntersection(pt a, pt b, pt c, pt d) { - if (orientation(a, b, c) == 0 && orientation(a, b, d) == 0) { - double dist = abs(a - b); - return (abs(a - c) <= dist && abs(b - c) <= dist) || - (abs(a - d) <= dist && abs(b - d) <= dist); - } - return orientation(a, b, c) * orientation(a, b, d) <= 0 && - orientation(c, d, a) * orientation(c, d, b) <= 0; -} - -// Berechnet die Schnittpunkte der Strecken a-b und c-d. Enthält entweder -// keinen Punkt, den einzigen Schnittpunkt oder die Endpunkte der -// Schnittstrecke. operator<, min, max müssen noch geschrieben werden! -vector<pt> lineSegmentIntersection(pt a, pt b, pt c, pt d) { - vector<pt> result; - if (orientation(a, b, c) == 0 && orientation(a, b, d) == 0 && - orientation(c, d, a) == 0 && orientation(c, d, b) == 0) { - pt minAB = min(a, b), maxAB = max(a, b); - pt minCD = min(c, d), maxCD = max(c, d); - if (minAB < minCD && maxAB < minCD) return result; - if (minCD < minAB && maxCD < minAB) return result; - pt start = max(minAB, minCD), end = min(maxAB, maxCD); - result.push_back(start); - if (start != end) result.push_back(end); - return result; - } - double x1 = real(b) - real(a), y1 = imag(b) - imag(a); - double x2 = real(d) - real(c), y2 = imag(d) - imag(c); - double u1 = (-y1 * (real(a) - real(c)) + x1 * (imag(a) - imag(c))) / - (-x2 * y1 + x1 * y2); - double u2 = (x2 * (imag(a) - imag(c)) - y2 * (real(a) - real(c))) / - (-x2 * y1 + x1 * y2); - if (u1 >= 0 && u1 <= 1 && u2 >= 0 && u2 <= 1) { - double x = real(a) + u2 * x1, y = imag(a) + u2 * y1; - result.push_back(pt(x, y)); - } - return result; -} - -// Entfernung von Punkt p zur Gearden durch a-b. -double distToLine(pt a, pt b, pt p) { - return abs(cross(p - a, b - a)) / abs(b - a); -} - -// Liegt p auf der Geraden a-b? -bool pointOnLine(pt a, pt b, pt p) { - return orientation(a, b, p) == 0; -} - -// Liegt p auf der Strecke a-b? -bool pointOnLineSegment(pt a, pt b, pt p) { - if (orientation(a, b, p) != 0) return false; - return real(p) >= min(real(a), real(b)) && - real(p) <= max(real(a), real(b)) && - imag(p) >= min(imag(a), imag(b)) && - imag(p) <= max(imag(a), imag(b)); -} - -// Entfernung von Punkt p zur Strecke a-b. -double distToSegment(pt a, pt b, pt p) { - if (a == b) return abs(p - a); - double segLength = abs(a - b); - double u = ((real(p) - real(a)) * (real(b) - real(a)) + - (imag(p) - imag(a)) * (imag(b) - imag(a))) / - (segLength * segLength); - pt projection(real(a) + u * (real(b) - real(a)), - imag(a) + u * (imag(b) - imag(a))); - double projectionDist = abs(p - projection); - if (!pointOnLineSegment(a, b, projection)) projectionDist = 1e30; - return min(projectionDist, min(abs(p - a), abs(p - b))); -} - -// Kürzeste Entfernung zwischen den Strecken a-b und c-d. -double distBetweenSegments(pt a, pt b, pt c, pt d) { - if (lineSegmentIntersection(a, b, c, d)) return 0.0; - double result = distToSegment(a, b, c); - result = min(result, distToSegment(a, b, d)); - result = min(result, distToSegment(c, d, a)); - return min(result, distToSegment(c, d, b)); + return (orien > EPS) - (orien < -EPS); } // Liegt d in der gleichen Ebene wie a, b, und c? bool isCoplanar(pt a, pt b, pt c, pt d) { - return abs((b - a) * (c - a) * (d - a)) < EPSILON; -} - -// Berechnet den Flächeninhalt eines Polygons (nicht selbstschneidend). -// Punkte gegen den Uhrzeigersinn: positiv, sonst negativ. -double areaOfPolygon(vector<pt> &polygon) { // Jeder Eckpunkt nur einmal. - double res = 0; int n = polygon.size(); - for (int i = 0; i < n; i++) - res += real(polygon[i]) * imag(polygon[(i + 1) % n]) - - real(polygon[(i + 1) % n]) * imag(polygon[i]); - return 0.5 * res; -} - -// Schneiden sich (p1, p2) und (p3, p4) (gegenüberliegende Ecken). -bool rectIntersection(pt p1, pt p2, pt p3, pt p4) { - double minx12=min(real(p1), real(p2)), maxx12=max(real(p1), real(p2)); - double minx34=min(real(p3), real(p4)), maxx34=max(real(p3), real(p4)); - double miny12=min(imag(p1), imag(p2)), maxy12=max(imag(p1), imag(p2)); - double miny34=min(imag(p3), imag(p4)), maxy34=max(imag(p3), imag(p4)); - return (maxx12 >= minx34) && (maxx34 >= minx12) && - (maxy12 >= miny34) && (maxy34 >= miny12); -} - -// Testet, ob ein Punkt im Polygon liegt (beliebige Polygone). -bool pointInPolygon(pt p, vector<pt> &polygon) { // Punkte nur einmal. - pt rayEnd = p + pt(1, 1000000); - int counter = 0, n = polygon.size(); - for (int i = 0; i < n; i++) { - pt start = polygon[i], end = polygon[(i + 1) % n]; - if (lineSegmentIntersection(p, rayEnd, start, end)) counter++; - } - return counter & 1; + return abs((b - a) * (c - a) * (d - a)) < EPS; } diff --git a/geometry/formulars3d.cpp b/geometry/formulars3d.cpp new file mode 100644 index 0000000..f527d57 --- /dev/null +++ b/geometry/formulars3d.cpp @@ -0,0 +1,53 @@ +// Skalarprodukt +double operator|(pt3 a, pt3 b) { + return a.x * b.x + a.y*b.y + a.z*b.z; +} +double dot(pt3 a, pt3 b) {return a|b;} + +// Kreuzprodukt +pt3 operator*(pt3 a, pt3 b) {return {a.y*b.z - a.z*b.y, + a.z*b.x - a.x*b.z, + a.x*b.y - a.y*b.x};} +pt3 cross(pt3 a, pt3 b) {return a*b;} + +// Länge von a +double abs(pt3 a) {return sqrt(dot(a, a));} +double abs(pt3 a, pt3 b) {return abs(b - a);} + +// Mixedprodukt +double mixed(pt3 a, pt3 b, pt3 c) {return a*b|c;}; + +// orientierung von p zu der Ebene durch a, b, c +// -1 => gegen den Uhrzeigersinn, +// 0 => kolliniear, +// 1 => im Uhrzeigersinn. +int orientation(pt3 a, pt3 b, pt3 c, pt3 p) { + double orien = mixed(b - a, c - a, p - a); + return (orien > EPS) - (orien < -EPS); +} + +// Entfernung von Punkt p zur Ebene a,b,c. +double distToPlane(pt3 a, pt3 b, pt3 c, pt3 p) { + pt3 n = cross(b-a, c-a); + return (abs(dot(n, p)) - dot(n, a)) / abs(n); +} + +// Liegt p in der Ebene a,b,c? +bool pointOnPlane(pt3 a, pt3 b, pt3 c, pt3 p) { + return orientation(a, b, c, p) == 0; +} + +// Schnittpunkt von der Grade a-b und der Ebene c,d,e +// die Grade darf nicht parallel zu der Ebene sein! +pt3 linePlaneIntersection(pt3 a, pt3 b, pt3 c, pt3 d, pt3 e) { + pt3 n = cross(d-c, e-c); + pt3 d = b - a; + return a - d * (dot(n, a) - dot(n, c)) / dot(n, d); +} + +// Abstand zwischen der Grade a-b und c-d +double lineLineDist(pt3 a, pt3 b, pt3 c, pt3 d) { + pt3 n = cross(b - a, d - c); + if (abs(n) < EPS) return distToLine(a, b, c); + return abs(dot(a - c, n)) / abs(n); +}
\ No newline at end of file diff --git a/geometry/geometry.tex b/geometry/geometry.tex index de7b2f6..1201917 100644 --- a/geometry/geometry.tex +++ b/geometry/geometry.tex @@ -1,13 +1,47 @@ \section{Geometrie} -\subsection{Closest Pair} -\lstinputlisting{geometry/closestPair.cpp} +\begin{algorithm}{Closest Pair} + \begin{methods} + \method{shortestDist}{kürzester Abstand zwischen Punkten}{n\*\log(n)} + \end{methods} + \sourcecode{geometry/closestPair.cpp} +\end{algorithm} -\subsection{Geraden} -\lstinputlisting{geometry/lines.cpp} +\begin{algorithm}{Konvexe Hülle} + \begin{methods} + \method{convexHull}{berechnet Konvexehülle}{n\*\log(n)} + \end{methods} + \begin{itemize} + \item Konvexehülle gegen den Uhrzeigersinn Sortiert + \item nur Eckpunkte enthalten(für alle Punkte = im CCW Test entfernen) + \item Erster und Letzter Punkt sind identisch + \end{itemize} + \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/triangle.cpp} +\sourcecode{geometry/polygon.cpp} +\sourcecode{geometry/circle.cpp} +\sourcecode{geometry/sortAround.cpp} -\subsection{Konvexe Hülle} -\lstinputlisting{geometry/convexHull.cpp} +\subsection{Formeln - 3D} +\sourcecode{geometry/formulars3d.cpp} -\subsection{Formeln - \lstinline{std::complex}} -\lstinputlisting{geometry/formulars.cpp} +\subsection{3D-Kugeln} +\sourcecode{geometry/spheres.cpp} + +\optional{ +\subsection{Geraden} +\sourcecode{geometry/lines.cpp} +} diff --git a/geometry/lines.cpp b/geometry/lines.cpp index 2594ab4..95536a4 100644 --- a/geometry/lines.cpp +++ b/geometry/lines.cpp @@ -1,32 +1,33 @@ -// Nicht complex<double> benutzen. Eigene struct schreiben. struct line { double a, b, c; // ax + by + c = 0; vertikale Line: b = 0, sonst: b = 1 + line(pt p, pt q) : a(-imag(q-p)), b(real(q-p)), c(cross({b, -a},p)) {} }; line pointsToLine(pt p1, pt p2) { line l; - if (fabs(p1.x - p2.x) < EPSILON) { - l.a = 1; l.b = 0.0; l.c = -p1.x; + if (abs(real(p1 - p2)) < EPS) { + l.a = 1; l.b = 0.0; l.c = -real(p1); } else { - l.a = -(double)(p1.y - p2.y) / (p1.x - p2.x); + l.a = -imag(p1 - p2) / real(p1 - p2); l.b = 1.0; - l.c = -(double)(l.a * p1.x) - p1.y; + l.c = -(l.a * real(p1)) - imag(p1); } return l; } -bool areParallel(line l1, line l2) { - return (fabs(l1.a - l2.a) < EPSILON) && (fabs(l1.b - l2.b) < EPSILON); +bool parallel(line l1, line l2) { + return (abs(l1.a - l2.a) < EPS) && (abs(l1.b - l2.b) < EPS); } -bool areSame(line l1, line l2) { - return areParallel(l1, l2) && (fabs(l1.c - l2.c) < EPSILON); +bool same(line l1, line l2) { + return parallel(l1, l2) && (abs(l1.c - l2.c) < EPS); } -bool areIntersect(line l1, line l2, pt &p) { - if (areParallel(l1, l2)) return false; - p.x = (l2.b * l1.c - l1.b * l2.c) / (l2.a * l1.b - l1.a * l2.b); - if (fabs(l1.b) > EPSILON) p.y = -(l1.a * p.x + l1.c); - else p.y = -(l2.a * p.x + l2.c); +bool intersect(line l1, line l2, pt& p) { + if (parallel(l1, l2)) return false; + double y, x = (l2.b * l1.c - l1.b * l2.c) / (l2.a * l1.b - l1.a * l2.b); + if (abs(l1.b) > EPS) y = -(l1.a * x + l1.c); + else y = -(l2.a * x + l2.c); + p = {x, y}; return true; } diff --git a/geometry/linesAndSegments.cpp b/geometry/linesAndSegments.cpp new file mode 100644 index 0000000..6c53cee --- /dev/null +++ b/geometry/linesAndSegments.cpp @@ -0,0 +1,90 @@ +// Test auf Streckenschnitt zwischen a-b und c-d. +bool lineSegmentIntersection(pt a, pt b, pt c, pt d) { + if (orientation(a, b, c) == 0 && orientation(a, b, d) == 0) + return pointOnLineSegment(a,b,c) || + pointOnLineSegment(a,b,d) || + pointOnLineSegment(c,d,a) || + pointOnLineSegment(c,d,b); + return orientation(a, b, c) * orientation(a, b, d) <= 0 && + orientation(c, d, a) * orientation(c, d, b) <= 0; +} + +// Berechnet die Schnittpunkte der Strecken p0-p1 und p2-p3. +// Enthält entweder keinen Punkt, den einzigen Schnittpunkt +// oder die Endpunkte der Schnittstrecke. +vector<pt> lineSegmentIntersection(pt p0, pt p1, pt p2, pt p3) { + double a = cross(p1 - p0, p3 - p2); + 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 (a > EPS) return {p0 + b/a*(p1 - p0)}; + vector<pt> result; + auto insertUnique = [&](pt p) { + for (auto q: result) if (abs(p - q) < EPS) return; + result.push_back(p); + }; + if (dot(p2-p0, p3-p0) < EPS) insertUnique(p0); + if (dot(p2-p1, p3-p1) < EPS) insertUnique(p1); + if (dot(p0-p2, p1-p2) < EPS) insertUnique(p2); + if (dot(p0-p3, p1-p3) < EPS) insertUnique(p3); + return result; +} + +// Entfernung von Punkt p zur Gearden durch a-b. 2d und 3d +double distToLine(pt a, pt b, pt p) { + return abs(cross(p - a, b - a)) / abs(b - a); +} + +// Liegt p auf der Geraden a-b? 2d und 3d +bool pointOnLine(pt a, pt b, pt p) { + return orientation(a, b, p) == 0; +} + +// Test auf Linienschnitt zwischen a-b und c-d. +bool lineIntersection(pt a, pt b, pt c, pt d) { + return abs(cross(a - b, c - d)) < EPS; +} + +// Berechnet den Schnittpunkt der Graden p0-p1 und p2-p3. +// die Graden dürfen nicht parallel sein! +pt lineIntersection(pt p0, pt p1, pt p2, pt p3) { + double a = cross(p1 - p0, p3 - p2); + double b = cross(p2 - p0, p3 - p2); + return {p0 + b/a*(p1 - p0)}; +} + +// Liegt p auf der Strecke a-b? +bool pointOnLineSegment(pt a, pt b, pt p) { + if (orientation(a, b, p) != 0) return false; + ld dist = norm(a - b); + return norm(a - p) <= dist && norm(b - p) <= dist; +} + +// Entfernung von Punkt p zur Strecke a-b. +double distToSegment(pt a, pt b, pt p) { + if (a == b) return abs(p - a); + pt dir = b - a; + if (dot(dir, a) <= dot(dir, p) && dot(dir, p) <= dot(dir, b)) { + return distToLine(a, b, p); + } else { + return min(abs(p - a), abs(p - b)); +}} + +// Kürzeste Entfernung zwischen den Strecken a-b und c-d. +double distBetweenSegments(pt a, pt b, pt c, pt d) { + if (lineSegmentIntersection(a, b, c, d)) return 0.0; + double result = distToSegment(a, b, c); + result = min(result, distToSegment(a, b, d)); + result = min(result, distToSegment(c, d, a)); + return min(result, distToSegment(c, d, b)); +} + +// sortiert alle Punkte pts auf einer Linie +// entsprechend der richtung dir 2d und 3d +void sortLine(pt dir, vector<pt>& pts) { + sort(pts.begin(), pts.end(), [&](pt a, pt b){ + return dot(dir, a) < dot(dir, b); + }); +}
\ No newline at end of file diff --git a/geometry/polygon.cpp b/geometry/polygon.cpp new file mode 100644 index 0000000..6a8080d --- /dev/null +++ b/geometry/polygon.cpp @@ -0,0 +1,33 @@ +// Berechnet den Flächeninhalt eines Polygons +// (nicht selbstschneidend). +// Punkte gegen den Uhrzeigersinn: positiv, sonst negativ. +double area(const vector<pt>& polygon) {//points unique + double res = 0; int n = sz(polygon); + for (int i = 0; i < n; i++) + res += cross(polygon[i], polygon[(i + 1) % n]); + return 0.5 * res; +} + +// Testet, ob ein Punkt im Polygon liegt (beliebige Polygone). +bool inside(pt p, const vector<pt>& polygon) {//points unique + pt rayEnd = p + pt(1, 1000000); + int counter = 0, n = sz(polygon); + for (int i = 0; i < n; i++) { + pt start = polygon[i], end = polygon[(i + 1) % n]; + if (lineSegmentIntersection(p, rayEnd, start, end)) { + counter++; + }} + return counter & 1; +} + +//convex hull without duplicates, h[0] == h.back() +bool inside(pt p, const vector<pt>& hull) { + if (orientation(hull[0], hull.back(), p) > 0) return false; + ll l = 0, r = sz(hull) - 1; + while (l + 1 < r) { + ll m = (l + r) / 2; + if (orientation(hull[0], hull[m], p) >= 0) l = m; + else r = m; + } + return orientation(hull[l], hull[r], p) >= 0; +} diff --git a/geometry/segmentIntersection.cpp b/geometry/segmentIntersection.cpp new file mode 100644 index 0000000..c2961da --- /dev/null +++ b/geometry/segmentIntersection.cpp @@ -0,0 +1,63 @@ +struct seg { + pt a, b; + int id; + bool operator<(const seg& o) const { + if (real(a) < real(o.a)) { + int s = orientation(a, b, o.a); + return (s > 0 || (s == 0 && imag(a) < imag(o.a)));// + } else if (real(a) > real(o.a)) { + int s = orientation(o.a, o.b, a); + return (s < 0 || (s == 0 && imag(a) < imag(o.a))); + } + return imag(a) < imag(o.a); + } +}; + +struct event { + pt p; + int id, type; + bool operator<(const event& o) const { + if (real(p) != real(o.p)) return real(p) < real(o.p); + if (type != o.type) return type > o.type; + return imag(p) < imag(o.p); + } +}; + +bool lessPT(const pt& a, const pt& b) { + return real(a) != real(b) ? real(a) < real(b) + : imag(a) < imag(b); +} + +bool intersect(const seg& a, const seg& b) { + return lineSegmentIntersection(a.a, a.b, b.a, b.b); +} + +pair<int, int> intersect(vector<seg>& segs) { + vector<event> events; + for (seg& s : segs) { + if (lessPT(s.b, s.a)) swap(s.b, s.a); + events.push_back({s.a, s.id, 1}); + events.push_back({s.b, s.id, -1}); + } + sort(all(events)); + + set<seg> q; + vector<set<seg>::iterator> where(sz(segs)); + for (auto e : events) { + int id = e.id; + if (e.type > 0) { + auto it = q.lower_bound(segs[id]); + if (it != q.end() && intersect(*it, segs[id])) + return {it->id, segs[id].id}; + if (it != q.begin() && intersect(*prev(it), segs[id])) + return {prev(it)->id, segs[id].id}; + where[id] = q.insert(it, segs[id]); + } else { + auto it = where[id]; + if (it != q.begin() && next(it) != q.end() && intersect(*next(it), *prev(it))) + return {next(it)->id, prev(it)->id}; + q.erase(it); + } + } + return {-1, -1}; +}
\ No newline at end of file diff --git a/geometry/sortAround.cpp b/geometry/sortAround.cpp new file mode 100644 index 0000000..a95d224 --- /dev/null +++ b/geometry/sortAround.cpp @@ -0,0 +1,10 @@ +bool left(pt p) {return real(p) < 0 ||
+ (real(p) == 0 && imag(p) < 0);}
+
+void sortAround(pt p, vector<pt>& ps) {
+ sort(all(ps), [&](const pt& a, const pt& b){
+ if (left(a - p) != left(b - p))
+ return left(a - p) > left(b - p);
+ return cross(p, a, b) > 0;
+ });
+}
\ No newline at end of file diff --git a/geometry/spheres.cpp b/geometry/spheres.cpp new file mode 100644 index 0000000..0657ab8 --- /dev/null +++ b/geometry/spheres.cpp @@ -0,0 +1,29 @@ +// Great Cirlce Distance mit Längen- und Breitengrad. +double gcDist(double pLat, double pLon, + double qLat, double qLon, double radius) { + pLat *= PI / 180; pLon *= PI / 180; + qLat *= PI / 180; qLon *= PI / 180; + return radius * acos(cos(pLat) * cos(pLon) * + cos(qLat) * cos(qLon) + + cos(pLat) * sin(pLon) * + cos(qLat) * sin(qLon) + + sin(pLat) * sin(qLat)); +} + +// Great Cirlce Distance mit kartesischen Koordinaten. +double gcDist(point p, point q) { + return acos(p.x * q.x + p.y * q.y + p.z * q.z); +} + +// 3D Punkt in kartesischen Koordinaten. +struct point{ + double x, y, z; + point() {} + point(double x, double y, double z) : x(x), y(y), z(z) {} + point(double lat, double lon) { + lat *= PI / 180.0; lon *= PI / 180.0; + x = cos(lat) * sin(lon); + y = cos(lat) * cos(lon); + z = sin(lat); + } +}; diff --git a/geometry/triangle.cpp b/geometry/triangle.cpp new file mode 100644 index 0000000..3a39302 --- /dev/null +++ b/geometry/triangle.cpp @@ -0,0 +1,40 @@ +// Mittelpunkt des Dreiecks abc. +pt centroid(pt a, pt b, pt c) {return (a + b + c) / 3.0;} + +// Flächeninhalt eines Dreicks bei bekannten Eckpunkten. +double area(pt a, pt b, pt c) { + return abs(cross(b - a, c - a)) / 2.0; +} + +// Flächeninhalt eines Dreiecks bei bekannten Seitenlängen. +double area(double a, double b, double c) { + double s = (a + b + c) / 2.0; + return sqrt(s * (s-a) * (s-b) * (s-c)); +} + +// Zentrum des Kreises durch alle Eckpunkte +pt outCenter(pt a, pt b, pt c) { + double d = 2.0 * (real(a) * imag(b-c) + + real(b) * imag(c-a) + + real(c) * imag(a-b)); + return (a*conj(a)*conj(b-c) + + b*conj(b)*conj(c-a) + + c*conj(c)*conj(a-b)) / d; +} + +// Zentrum des größten Kreises im Dreiecke +pt inCenter(pt a, pt b, pt c) { + double x = abs(a-b), y = abs(b-c), z = abs(a-c); + return (y*a + z*b + x*c) / (a+b+c); +} + + +// 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) + ); +} |