Skip to content

Add cross, dot and norm operators #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
bcs5 wants to merge 7 commits into
base: master
Choose a base branch
from
Draft

Add cross, dot and norm operators #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f718893
Add cross, dot and norm operators
bcs5 Aug 12, 2023
41cdc9a
Remove dist funcs
bcs5 Aug 12, 2023
8cbb5ba
nits
bcs5 Aug 12, 2023
d793eb2
Use operators on convex hull
bcs5 Aug 12, 2023
4cebf75
Fix projPtLine
bcs5 Aug 12, 2023
540d4d8
fixup! Use operators on convex hull
bcs5 Aug 13, 2023
27ae542
fixup! fixup! Use operators on convex hull
bcs5 Aug 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 8 additions & 8 deletions code/Geometry/ConvexHull.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,11 @@ vector<PT> convexHull(vector<PT> p, bool needs = 1) {
if(n <= 1) return p;
vector<PT> 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
Expand Down Expand Up @@ -40,28 +40,28 @@ int maximizeScalarProduct(const vector<PT> &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;
}
}
Expand Down
110 changes: 51 additions & 59 deletions code/Geometry/Geometry.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,16 +1,23 @@
#include <bits/stdc++.h>
using namespace std;

const double inf = 1e100, eps = 1e-9;
const double PI = acos(-1.0L);
int cmp (double a, double b = 0) {
if (abs(a-b) < eps) return 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;
Expand All @@ -21,13 +28,8 @@ 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 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);
Expand All @@ -39,69 +41,57 @@ 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<PT, int> 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);
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 = dot(b-a, b-a);
b = b - a, c = c - a;
double r = b * b;
if (cmp(r) == 0) return a;
r = dot(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
return dist(c, projPtSeg(a, b, c));
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
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) {
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) {
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
// Para intersecao de segmentos, cheque primeiro se os segmentos se intersectam e que nao sao paralelos
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
Expand All @@ -110,7 +100,7 @@ PT circleCenter (PT a, PT b, PT c) {
}
vector<PT> circle2PtsRad (PT p1, PT p2, double r) {
vector<PT> 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);
Expand All @@ -122,48 +112,48 @@ vector<PT> 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;
bool circleLineIntersection (PT a, PT b, PT c, double r) {
return cmp(!(c - projPtLine(a, b, c)), r) <= 0;
}
vector<PT> circleLine (PT a, PT b, PT c, double r) {
vector<PT> 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);
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;
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<PT> &p, PT q) {
bool pointInConvexPolygon (const vector<PT> &p, PT q) {
if (p.size() == 1) return p.front() == 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]);
}
// 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<PT> &p, PT q) {
bool pointInPolygon (const vector<PT> &p, PT q) {
bool c = 0;
for(int i = 0; i < p.size(); i++){
int j = (i + 1) % p.size();
Expand All @@ -175,12 +165,12 @@ bool pointInPolygon(const vector<PT> &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<PT> circleCircle (PT a, double r, PT b, double R) {
vector<PT> 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);
Expand All @@ -204,8 +194,8 @@ double computeSignedArea (const vector<PT> &p) {
}
return area/2.0;
}
double computeArea(const vector<PT> &p) { return abs(computeSignedArea(p)); }
PT computeCentroid(const vector<PT> &p) {
double computeArea (const vector<PT> &p) { return abs(computeSignedArea(p)); }
PT computeCentroid (const vector<PT> &p) {
PT c(0,0);
double scale = 6.0 * computeSignedArea(p);
for(int i = 0; i < p.size(); i++){
Expand All @@ -215,7 +205,7 @@ PT computeCentroid(const vector<PT> &p) {
return c / scale;
}
// Testa se o poligno listada em ordem CW ou CCW eh simples (nenhuma linha se intersecta)
bool isSimple(const vector<PT> &p) {
bool isSimple (const vector<PT> &p) {
for(int i = 0; i < p.size(); i++) {
for(int k = i + 1; k < p.size(); k++) {
int j = (i + 1) % p.size();
Expand All @@ -227,10 +217,12 @@ bool isSimple(const vector<PT> &p) {
}
return true;
}

PT normalize (PT p) { return p/!p; }
vector< pair<PT, PT> > getTangentSegs (PT c1, double r1, PT c2, double r2) {
if (r1 < r2) swap(c1, c2), swap(r1, r2);
vector<pair<PT, PT> > 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;
Expand Down