From f71889365f80dea63fc540a97c5c9bf85f0b0df4 Mon Sep 17 00:00:00 2001 From: Bezaliel Silva Date: Sat, 12 Aug 2023 09:14:29 -0300 Subject: [PATCH 1/7] Add cross, dot and norm operators --- code/Geometry/Geometry.cpp | 78 +++++++++++++++++--------------------- 1 file changed, 35 insertions(+), 43 deletions(-) diff --git a/code/Geometry/Geometry.cpp b/code/Geometry/Geometry.cpp index 312ae3d..2d08671 100644 --- a/code/Geometry/Geometry.cpp +++ b/code/Geometry/Geometry.cpp @@ -1,3 +1,6 @@ +#include +using namespace std; + const double inf = 1e100, eps = 1e-9; const double PI = acos(-1.0L); int cmp (double a, double b = 0) { @@ -5,12 +8,16 @@ int cmp (double a, double b = 0) { return (a < b) ? -1 : +1; } struct PT { - double x, y; - PT(double x = 0, double y = 0) : x(x), y(y) {} + typedef long long T; + T x, y; + PT(T x = 0, T y = 0) : x(x), y(y) {} PT operator + (const PT &p) const { return PT(x+p.x, y+p.y); } PT operator - (const PT &p) const { return PT(x-p.x, y-p.y); } PT operator * (double c) const { return PT(x*c, y*c); } PT operator / (double c) const { return PT(x/c, y/c); } + T operator * (const PT &p) const { return x*p.x + y*p.y; } + T operator % (const PT &p) const { return x*p.y - y*p.x; } + double operator !() const { return hypot(x, y); } bool operator <(const PT &p) const { if(cmp(x, p.x) != 0) return x < p.x; return cmp(y, p.y) < 0; @@ -21,13 +28,10 @@ struct PT { ostream &operator<<(ostream &os, const PT &p) { return os << "(" << p.x << "," << p.y << ")"; } -double dot (PT p, PT q) { return p.x * q.x + p.y*q.y; } -double cross (PT p, PT q) { return p.x * q.y - p.y*q.x; } -double dist2 (PT p, PT q = PT(0, 0)) { return dot(p-q, p-q); } + +double dist2 (PT p, PT q = PT(0, 0)) { return (p-q)*(p-q); } double dist (PT p, PT q) { return hypot(p.x-q.x, p.y-q.y); } -double norm (PT p) { return hypot(p.x, p.y); } -PT normalize (PT p) { return p/hypot(p.x, p.y); } -double angle (PT p, PT q) { return atan2(cross(p, q), dot(p, q)); } +double angle (PT p, PT q) { return atan2(p % q, p * q); } double angle (PT p) { return atan2(p.y, p.x); } double polarAngle (PT p) { double a = atan2(p.y,p.x); @@ -39,31 +43,17 @@ PT rotateCW90 (PT p) { return PT(p.y, -p.x); } PT rotateCCW (PT p, double t) { return PT(p.x*cos(t)-p.y*sin(t), p.x*sin(t)+p.y*cos(t)); } -typedef pair Line; -PT getDir (PT a, PT b) { - if (a.x == b.x) return PT(0, 1); - if (a.y == b.y) return PT(1, 0); - int dx = b.x-a.x; - int dy = b.y-a.y; - int g = __gcd(abs(dx), abs(dy)); - if (dx < 0) g = -g; - return PT(dx/g, dy/g); -} -Line getLine (PT a, PT b) { - PT dir = getDir(a, b); - return {dir, cross(dir, a)}; -} PT projPtLine (PT a, PT b, PT c) { // ponto c na linha a - b, a.b = |a| cost * |b| - return a + (b-a) * dot(b-a, c-a)/dot(b-a, b-a); + return a + (b-a) * ((b-a)*(c-a))/((b-a)*(b-a)); } PT reflectPointLine (PT a, PT b, PT c) { PT p = projPtLine(a, b, c); return p*2 - c; } PT projPtSeg (PT a, PT b, PT c) { // c no segmento a - b - double r = dot(b-a, b-a); + double r = (b-a)*(b-a); if (cmp(r) == 0) return a; - r = dot(b-a, c-a)/r; + r = (b-a)*(c-a)/r; if (cmp(r, 0) < 0) return a; if (cmp(r, 1) > 0) return b; return a + (b - a) * r; @@ -74,13 +64,13 @@ double distancePointSegment (PT a, PT b, PT c) { // ponto c e o segmento a - b bool ptInSegment (PT a, PT b, PT c) { // ponto c esta em um segmento a - b if (a == b) return a == c; a = a-c, b = b-c; - return cmp(cross(a, b)) == 0 && cmp(dot(a, b)) <= 0; + return cmp(a % b) == 0 && cmp(a * b) <= 0; } bool parallel (PT a, PT b, PT c, PT d) { - return cmp(cross(b - a, c - d)) == 0; + return cmp((b - a)%(c - d)) == 0; } bool collinear (PT a, PT b, PT c, PT d) { - return parallel(a, b, c, d) && cmp(cross(a - b, a - c)) == 0 && cmp(cross(c - d, c - a)) == 0; + return parallel(a, b, c, d) && cmp((a - b)%(a - c)) == 0 && cmp((c - d)%(c - a)) == 0; } // Calcula distancia entre o ponto (x, y, z) e o plano ax + by + cz = d double distPtPlane(double x, double y, double z, double a, double b, double c, double d) { @@ -89,11 +79,11 @@ double distPtPlane(double x, double y, double z, double a, double b, double c, d bool segInter (PT a, PT b, PT c, PT d) { if (collinear(a, b, c, d)) { if (a == c || a == d || b == c || b == d) return true; - if (cmp(dot(c - a, c - b)) > 0 && cmp(dot(d - a, d - b)) > 0 && cmp(dot(c - b, d - b)) > 0) return false; + if (cmp((c - a)*(c - b)) > 0 && cmp((d - a)*(d - b)) > 0 && cmp((c - b)*(d - b)) > 0) return false; return true; } - if (cmp(cross(d - a, b - a) * cross(c - a, b - a)) > 0) return false; - if (cmp(cross(a - c, d - c) * cross(b - c, d - c)) > 0) return false; + if (cmp((d - a)%(b - a) * ((c - a)%(b - a))) > 0) return false; + if (cmp((a - c)%(d - c) * ((b - c)%(d - c))) > 0) return false; return true; } // Calcula a intersecao entre as retas a - b e c - d assumindo que uma unica intersecao existe @@ -101,7 +91,7 @@ bool segInter (PT a, PT b, PT c, PT d) { PT lineLine (PT a, PT b, PT c, PT d) { b = b - a; d = c - d; c = c - a; // assert(cmp(cross(b, d)) != 0); - return a + b * cross(c, d) / cross(b, d); + return a + b * (c % d) / (b % d); } PT circleCenter (PT a, PT b, PT c) { b = (a + b) / 2; // bissector @@ -128,26 +118,26 @@ bool circleLineIntersection(PT a, PT b, PT c, double r) { vector circleLine (PT a, PT b, PT c, double r) { vector ret; PT p = projPtLine(a, b, c), p1; - double h = norm(c-p); + double h = !(c-p); if (cmp(h,r) == 0) { ret.push_back(p); } else if (cmp(h,r) < 0) { double k = sqrt(r*r - h*h); - p1 = p + (b-a)/(norm(b-a))*k; + p1 = p + (b-a)/(!(b-a))*k; ret.push_back(p1); - p1 = p - (b-a)/(norm(b-a))*k; + p1 = p - (b-a)/(!(b-a))*k; ret.push_back(p1); } return ret; } bool ptInsideTriangle(PT p, PT a, PT b, PT c) { - if(cross(b-a, c-b) < 0) swap(a, b); + if((b-a) % (c-b) < 0) swap(a, b); if(ptInSegment(a,b,p)) return 1; if(ptInSegment(b,c,p)) return 1; if(ptInSegment(c,a,p)) return 1; - bool x = cross(b-a, p-b) < 0; - bool y = cross(c-b, p-c) < 0; - bool z = cross(a-c, p-a) < 0; + bool x = (b-a)%(p-b) < 0; + bool y = (c-b)%(p-c) < 0; + bool z = (a-c)%(p-a) < 0; return x == y && y == z; } bool pointInConvexPolygon(const vector &p, PT q) { @@ -155,7 +145,7 @@ bool pointInConvexPolygon(const vector &p, PT q) { int l = 1, r = p.size()-1; while(abs(r-l) > 1) { int m = (r+l)/2; - if(cross(p[m]-p[0] , q-p[0]) < 0) r = m; + if((p[m]-p[0])%(q-p[0]) < 0) r = m; else l = m; } return ptInsideTriangle(q, p[0], p[l], p[r]); @@ -175,12 +165,12 @@ bool pointInPolygon(const vector &p, PT q) { } // area / semiperimeter double rIncircle (PT a, PT b, PT c) { - double ab = norm(a-b), bc = norm(b-c), ca = norm(c-a); - return abs(cross(b-a, c-a)/(ab+bc+ca)); + double ab = !(a-b), bc = !(b-c), ca = !(c-a); + return abs((b-a)%(c-a)/(ab+bc+ca)); } vector circleCircle (PT a, double r, PT b, double R) { vector ret; - double d = norm(a-b); + double d = !(a-b); if (d > r + R || d + min(r, R) < max(r, R)) return ret; double x = (d*d - R*R + r*r) / (2*d); // x = r*cos(R opposite angle) double y = sqrt(r*r - x*x); @@ -227,6 +217,8 @@ bool isSimple(const vector &p) { } return true; } + +PT normalize (PT p) { return p/!p; } vector< pair > getTangentSegs (PT c1, double r1, PT c2, double r2) { if (r1 < r2) swap(c1, c2), swap(r1, r2); vector > ans; From 41cdc9a8adebd03e4728deaf7d0acdbe76cb770a Mon Sep 17 00:00:00 2001 From: Bezaliel Silva Date: Sat, 12 Aug 2023 09:58:03 -0300 Subject: [PATCH 2/7] Remove dist funcs --- code/Geometry/Geometry.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/code/Geometry/Geometry.cpp b/code/Geometry/Geometry.cpp index 2d08671..7e65ea9 100644 --- a/code/Geometry/Geometry.cpp +++ b/code/Geometry/Geometry.cpp @@ -29,8 +29,6 @@ ostream &operator<<(ostream &os, const PT &p) { return os << "(" << p.x << "," << p.y << ")"; } -double dist2 (PT p, PT q = PT(0, 0)) { return (p-q)*(p-q); } -double dist (PT p, PT q) { return hypot(p.x-q.x, p.y-q.y); } double angle (PT p, PT q) { return atan2(p % q, p * q); } double angle (PT p) { return atan2(p.y, p.x); } double polarAngle (PT p) { @@ -59,7 +57,7 @@ PT projPtSeg (PT a, PT b, PT c) { // c no segmento a - b return a + (b - a) * r; } double distancePointSegment (PT a, PT b, PT c) { // ponto c e o segmento a - b - return dist(c, projPtSeg(a, b, c)); + return !(c - projPtSeg(a, b, c)); } bool ptInSegment (PT a, PT b, PT c) { // ponto c esta em um segmento a - b if (a == b) return a == c; @@ -100,7 +98,7 @@ PT circleCenter (PT a, PT b, PT c) { } vector circle2PtsRad (PT p1, PT p2, double r) { vector ret; - double d2 = dist2(p1, p2); + double d2 = p1 * p2; double det = r * r / d2 - 0.25; if (det < 0.0) return ret; double h = sqrt(det); @@ -113,7 +111,7 @@ vector circle2PtsRad (PT p1, PT p2, double r) { return ret; } bool circleLineIntersection(PT a, PT b, PT c, double r) { - return cmp(dist(c, projPtLine(a, b, c)), r) <= 0; + return cmp(!(c - projPtLine(a, b, c)), r) <= 0; } vector circleLine (PT a, PT b, PT c, double r) { vector ret; @@ -222,7 +220,7 @@ PT normalize (PT p) { return p/!p; } vector< pair > getTangentSegs (PT c1, double r1, PT c2, double r2) { if (r1 < r2) swap(c1, c2), swap(r1, r2); vector > ans; - double d = dist(c1, c2); + double d = !(c1 - c2); if (cmp(d) <= 0) return ans; double dr = abs(r1 - r2), sr = r1 + r2; if (cmp(dr, d) >= 0) return ans; From 8cbb5baf13cb7f84c86ca73729409d8e9b93a853 Mon Sep 17 00:00:00 2001 From: Bezaliel Silva Date: Sat, 12 Aug 2023 12:05:12 -0300 Subject: [PATCH 3/7] nits --- code/Geometry/Geometry.cpp | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/code/Geometry/Geometry.cpp b/code/Geometry/Geometry.cpp index 7e65ea9..51844ad 100644 --- a/code/Geometry/Geometry.cpp +++ b/code/Geometry/Geometry.cpp @@ -42,21 +42,23 @@ PT rotateCCW (PT p, double t) { return PT(p.x*cos(t)-p.y*sin(t), p.x*sin(t)+p.y*cos(t)); } PT projPtLine (PT a, PT b, PT c) { // ponto c na linha a - b, a.b = |a| cost * |b| - return a + (b-a) * ((b-a)*(c-a))/((b-a)*(b-a)); + PT b = b - a, c = c - a; + return a + b*(c*b)/(b*b); } PT reflectPointLine (PT a, PT b, PT c) { PT p = projPtLine(a, b, c); return p*2 - c; } PT projPtSeg (PT a, PT b, PT c) { // c no segmento a - b - double r = (b-a)*(b-a); + b = b - a, c = c - a; + double r = b * b; if (cmp(r) == 0) return a; - r = (b-a)*(c-a)/r; + r = c * b / r; if (cmp(r, 0) < 0) return a; - if (cmp(r, 1) > 0) return b; - return a + (b - a) * r; + if (cmp(r, 1) > 0) return a + b; + return a + b * r; } -double distancePointSegment (PT a, PT b, PT c) { // ponto c e o segmento a - b +double distPtSeg (PT a, PT b, PT c) { // ponto c e o segmento a - b return !(c - projPtSeg(a, b, c)); } bool ptInSegment (PT a, PT b, PT c) { // ponto c esta em um segmento a - b @@ -71,7 +73,7 @@ bool collinear (PT a, PT b, PT c, PT d) { return parallel(a, b, c, d) && cmp((a - b)%(a - c)) == 0 && cmp((c - d)%(c - a)) == 0; } // Calcula distancia entre o ponto (x, y, z) e o plano ax + by + cz = d -double distPtPlane(double x, double y, double z, double a, double b, double c, double d) { +double distPtPlane (double x, double y, double z, double a, double b, double c, double d) { return abs(a * x + b * y + c * z - d) / sqrt(a * a + b * b + c * c); } bool segInter (PT a, PT b, PT c, PT d) { @@ -110,7 +112,7 @@ vector circle2PtsRad (PT p1, PT p2, double r) { } return ret; } -bool circleLineIntersection(PT a, PT b, PT c, double r) { +bool circleLineIntersection (PT a, PT b, PT c, double r) { return cmp(!(c - projPtLine(a, b, c)), r) <= 0; } vector circleLine (PT a, PT b, PT c, double r) { @@ -128,7 +130,7 @@ vector circleLine (PT a, PT b, PT c, double r) { } return ret; } -bool ptInsideTriangle(PT p, PT a, PT b, PT c) { +bool ptInsideTriangle (PT p, PT a, PT b, PT c) { if((b-a) % (c-b) < 0) swap(a, b); if(ptInSegment(a,b,p)) return 1; if(ptInSegment(b,c,p)) return 1; @@ -138,7 +140,7 @@ bool ptInsideTriangle(PT p, PT a, PT b, PT c) { bool z = (a-c)%(p-a) < 0; return x == y && y == z; } -bool pointInConvexPolygon(const vector &p, PT q) { +bool pointInConvexPolygon (const vector &p, PT q) { if (p.size() == 1) return p.front() == q; int l = 1, r = p.size()-1; while(abs(r-l) > 1) { @@ -151,7 +153,7 @@ bool pointInConvexPolygon(const vector &p, PT q) { // Determina se o ponto esta num poligono possivelmente nao-convexo // Retorna 1 para pontos estritamente dentro, 0 para pontos estritamente fora do poligno // e 0 ou 1 para os pontos restantes -bool pointInPolygon(const vector &p, PT q) { +bool pointInPolygon (const vector &p, PT q) { bool c = 0; for(int i = 0; i < p.size(); i++){ int j = (i + 1) % p.size(); @@ -192,8 +194,8 @@ double computeSignedArea (const vector &p) { } return area/2.0; } -double computeArea(const vector &p) { return abs(computeSignedArea(p)); } -PT computeCentroid(const vector &p) { +double computeArea (const vector &p) { return abs(computeSignedArea(p)); } +PT computeCentroid (const vector &p) { PT c(0,0); double scale = 6.0 * computeSignedArea(p); for(int i = 0; i < p.size(); i++){ @@ -203,7 +205,7 @@ PT computeCentroid(const vector &p) { return c / scale; } // Testa se o poligno listada em ordem CW ou CCW eh simples (nenhuma linha se intersecta) -bool isSimple(const vector &p) { +bool isSimple (const vector &p) { for(int i = 0; i < p.size(); i++) { for(int k = i + 1; k < p.size(); k++) { int j = (i + 1) % p.size(); From d793eb20cd1bcdb94f92ac3105f1fac2fcc9713a Mon Sep 17 00:00:00 2001 From: Bezaliel Silva Date: Sat, 12 Aug 2023 20:55:31 -0300 Subject: [PATCH 4/7] Use operators on convex hull --- code/Geometry/ConvexHull.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/code/Geometry/ConvexHull.cpp b/code/Geometry/ConvexHull.cpp index 3e5a27b..e2f0c77 100644 --- a/code/Geometry/ConvexHull.cpp +++ b/code/Geometry/ConvexHull.cpp @@ -5,11 +5,11 @@ vector convexHull(vector p, bool needs = 1) { if(n <= 1) return p; vector h(n + 2); // se der wa bota n*2 for(int i = 0; i < n; i++) { - while(k >= 2 && cross(h[k - 1] - h[k - 2], p[i] - h[k - 2]) <= 0) k--; + while(k >= 2 && (h[k - 1] - h[k - 2]) % (p[i] - h[k - 2]) <= 0) k--; h[k++] = p[i]; } for(int i = n - 2, t = k + 1; i >= 0; i--) { - while(k >= t && cross(h[k - 1] - h[k - 2], p[i] - h[k - 2]) <= 0) k--; + while(k >= t && (h[k - 1] - h[k - 2]) % (p[i] - h[k - 2]) <= 0) k--; h[k++] = p[i]; } h.resize(k); // n+1 points where the first is equal to the last @@ -40,28 +40,28 @@ int maximizeScalarProduct(const vector &hull, PT vec) { int n = hull.size(); if(n < 20) { for(int i = 0; i < n; i++) { - if(dot(hull[i], vec) > dot(hull[ans], vec)) { + if(hull[i] * vec > hull[ans] * vec) { ans = i; } } } else { - if(dot(hull[1], vec) > dot(hull[ans], vec)) { + if(hull[1] * vec > hull[ans] * vec)) { ans = 1; } for(int rep = 0; rep < 2; rep++) { int l = 2, r = n - 1; while(l != r) { int mid = (l + r + 1) / 2; - bool flag = dot(hull[mid], vec) >= dot(hull[mid-1], vec); - if(rep == 0) { flag = flag && dot(hull[mid], vec) >= dot(hull[0], vec); } - else { flag = flag || dot(hull[mid-1], vec) < dot(hull[0], vec); } + bool flag = hull[mid] * vec >= hull[mid-1] * vec; + if(rep == 0) { flag = flag && hull[mid] * vec >= hull[0] * vec; } + else { flag = flag || hull[mid-1] * vec < hull[0] * vec; } if(flag) { l = mid; } else { r = mid - 1; } } - if(dot(hull[ans], vec) < dot(hull[l], vec)) { + if(hull[ans] * vec) < hull[l] * vec) { ans = l; } } From 4cebf75a6c8f5a5e6e229ae506768e8d25e54d70 Mon Sep 17 00:00:00 2001 From: Bezaliel Silva Date: Sat, 12 Aug 2023 20:58:15 -0300 Subject: [PATCH 5/7] Fix projPtLine --- code/Geometry/Geometry.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/Geometry/Geometry.cpp b/code/Geometry/Geometry.cpp index 51844ad..56d3a15 100644 --- a/code/Geometry/Geometry.cpp +++ b/code/Geometry/Geometry.cpp @@ -42,7 +42,7 @@ PT rotateCCW (PT p, double t) { return PT(p.x*cos(t)-p.y*sin(t), p.x*sin(t)+p.y*cos(t)); } PT projPtLine (PT a, PT b, PT c) { // ponto c na linha a - b, a.b = |a| cost * |b| - PT b = b - a, c = c - a; + b = b - a, c = c - a; return a + b*(c*b)/(b*b); } PT reflectPointLine (PT a, PT b, PT c) { From 540d4d8d06f308a5cc45d2c9dfcf50f3d4ae23cc Mon Sep 17 00:00:00 2001 From: Bezaliel Silva Date: Sat, 12 Aug 2023 22:03:48 -0300 Subject: [PATCH 6/7] fixup! Use operators on convex hull --- code/Geometry/ConvexHull.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/Geometry/ConvexHull.cpp b/code/Geometry/ConvexHull.cpp index e2f0c77..15ffc55 100644 --- a/code/Geometry/ConvexHull.cpp +++ b/code/Geometry/ConvexHull.cpp @@ -61,7 +61,7 @@ int maximizeScalarProduct(const vector &hull, PT vec) { r = mid - 1; } } - if(hull[ans] * vec) < hull[l] * vec) { + if(hull[ans] * vec < hull[l] * vec) { ans = l; } } From 27ae542723416bea014228ad4e0b851b7752a32d Mon Sep 17 00:00:00 2001 From: Bezaliel Silva Date: Sat, 12 Aug 2023 22:04:32 -0300 Subject: [PATCH 7/7] fixup! fixup! Use operators on convex hull --- code/Geometry/ConvexHull.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/Geometry/ConvexHull.cpp b/code/Geometry/ConvexHull.cpp index 15ffc55..5638e8f 100644 --- a/code/Geometry/ConvexHull.cpp +++ b/code/Geometry/ConvexHull.cpp @@ -45,7 +45,7 @@ int maximizeScalarProduct(const vector &hull, PT vec) { } } } else { - if(hull[1] * vec > hull[ans] * vec)) { + if(hull[1] * vec > hull[ans] * vec) { ans = 1; } for(int rep = 0; rep < 2; rep++) {