Skip to content

Commit

Permalink
geometry template update
Browse files Browse the repository at this point in the history
  • Loading branch information
Asad committed Mar 22, 2023
1 parent 007d6b6 commit 14ea361
Show file tree
Hide file tree
Showing 4 changed files with 306 additions and 162 deletions.
126 changes: 51 additions & 75 deletions geometry/convex-hull.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,78 +1,54 @@
// complexity O(n*logn)
// input: set of 2D points
// output: convex hull with sorted counter-clockwise

bool cmp(point a, point b) { return a.x < b.x || (a.x == b.x && a.y < b.y); }
bool cw(point a, point b, point c) { return a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y) < 0; }
bool ccw(point a, point b, point c) { return a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y) > 0; }

void convex_hull(int n, vector<point>& v) {
if (n == 1) return;
sort(all(v), &cmp);
if (v[0].x == v[n - 1].x && v[0].y == v[n - 1].y) {
point a = v[0];
v.clear();
v.PB(a);
return;
struct pt {
double x, y;
};

int orientation(pt a, pt b, pt c) {
double v = a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y);
if (v < 0) return -1; // clockwise
if (v > 0) return +1; // counter-clockwise
return 0;
}

bool cw(pt a, pt b, pt c, bool include_collinear) {
int o = orientation(a, b, c);
return o < 0 || (include_collinear && o == 0);
}
bool ccw(pt a, pt b, pt c, bool include_collinear) {
int o = orientation(a, b, c);
return o > 0 || (include_collinear && o == 0);
}

void convex_hull(vector<pt>& a, bool include_collinear = false) {
if (a.size() == 1)
return;

sort(a.begin(), a.end(), [](pt a, pt b) {
return make_pair(a.x, a.y) < make_pair(b.x, b.y);
});
pt p1 = a[0], p2 = a.back();
vector<pt> up, down;
up.push_back(p1);
down.push_back(p1);
for (int i = 1; i < (int)a.size(); i++) {
if (i == a.size() - 1 || cw(p1, a[i], p2, include_collinear)) {
while (up.size() >= 2 && !cw(up[up.size()-2], up[up.size()-1], a[i], include_collinear))
up.pop_back();
up.push_back(a[i]);
}
if (i == a.size() - 1 || ccw(p1, a[i], p2, include_collinear)) {
while (down.size() >= 2 && !ccw(down[down.size()-2], down[down.size()-1], a[i], include_collinear))
down.pop_back();
down.push_back(a[i]);
}
}
point p1 = v[0], p2 = v[n - 1];
vector<point> up, down;
up.PB(p1);
down.PB(p1);
for (int K = 1; K < n; K++) {
if (K == n - 1 || cw(p1, v[K], p2)) {
// while(sz(up)>=2 && !cw(up[sz(up)-2], up[sz(up)-1], v[K]))
// up.pop_back();
up.PB(v[K]);
}
if (K == n - 1 || ccw(p1, v[K], p2)) {
// while(sz(down)>=2 && !ccw(down[sz(down)-2],down[sz(down)-1], v[K]))
// down.pop_back();
down.PB(v[K]);
}

if (include_collinear && up.size() == a.size()) {
reverse(a.begin(), a.end());
return;
}
v.clear();
for (int K = 0; K < sz(down); K++) v.PB(down[K]);
for (int K = sz(up) - 2; K > 0; K--) v.PB(up[K]);
}
// from CP3
point pivot(0, 0);
bool angleCmp(point a, point b) { // angle-sorting function
if (collinear(pivot, a, b)) // special case
return dist(pivot, a) < dist(pivot, b); // check which one is closer
double d1x = a.x - pivot.x, d1y = a.y - pivot.y;
double d2x = b.x - pivot.x, d2y = b.y - pivot.y;
return (atan2(d1y, d1x) - atan2(d2y, d2x)) < 0;
} // compare two angles
vector<point> CH(vector<point> P) { // the content of P may be reshuffled
int i, j, n = (int)P.size();
if (n <= 3) {
if (!(P[0] == P[n - 1])) P.push_back(P[0]); // safeguard from corner case
return P;
} // special case, the CH is P itself
// first, find P0 = point with lowest Y and if tie: rightmost X
int P0 = 0;
for (i = 1; i < n; i++)
if (P[i].y < P[P0].y || (P[i].y == P[P0].y && P[i].x > P[P0].x)) P0 = i;
point temp = P[0];
P[0] = P[P0];
P[P0] = temp; // swap P[P0] with P[0]
// second, sort points by angle w.r.t. pivot P0
pivot = P[0]; // use this global variable as reference
sort(++P.begin(), P.end(), angleCmp); // we do not sort P[0]
// third, the ccw tests
vector<point> S;
S.push_back(P[n - 1]);
S.push_back(P[0]);
S.push_back(P[1]); // initial S
i = 2; // then, we check the rest
while (i < n) { // note: N must be >= 3 for this method to work
j = (int)S.size() - 1;
if (ccw(S[j - 1], S[j], P[i]))
S.push_back(P[i++]); // left turn, accept
else
S.pop_back();
} // or pop the top of S until we have a left turn
return S;
a.clear();
for (int i = 0; i < (int)up.size(); i++)
a.push_back(up[i]);
for (int i = down.size() - 2; i > 0; i--)
a.push_back(down[i]);
}

195 changes: 195 additions & 0 deletions geometry/line.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,195 @@
struct line {
PT a, b; // goes through points a and b
PT v; double c; //line form: direction vec [cross] (x, y) = c
line() {}
//direction vector v and offset c
line(PT v, double c) : v(v), c(c) {
auto p = get_points();
a = p.first; b = p.second;
}
// equation ax + by + c = 0
line(double _a, double _b, double _c) : v({_b, -_a}), c(-_c) {
auto p = get_points();
a = p.first; b = p.second;
}
// goes through points p and q
line(PT p, PT q) : v(q - p), c(cross(v, p)), a(p), b(q) {}
pair<PT, PT> get_points() { //extract any two points from this line
PT p, q; double a = -v.y, b = v.x; // ax + by = -c
if (sign(a) == 0) {
p = PT(0, -c / b);
q = PT(1, -c / b);
}
else if (sign(b) == 0) {
p = PT(-c / a, 0);
q = PT(-c / a, 1);
}
else {
p = PT(0, -c / b);
q = PT(1, (-c - a) / b);
}
return {p, q};
}
//ax + by + c = 0
array<double, 3> get_abc() {
double a = -v.y, b = v.x;
return {a, b, c};
}
// 1 if on the left, -1 if on the right, 0 if on the line
int side(PT p) { return sign(cross(v, p) - c); }
// line that is perpendicular to this and goes through point p
line perpendicular_through(PT p) { return {p, p + perp(v)}; }
// translate the line by vector t i.e. shifting it by vector t
line translate(PT t) { return {v, c + cross(v, t)}; }
// compare two points by their orthogonal projection on this line
// a projection point comes before another if it comes first according to vector v
bool cmp_by_projection(PT p, PT q) { return dot(v, p) < dot(v, q); }
line shift_left(double d) {
PT z = v.perp().truncate(d);
return line(a + z, b + z);
}
};
// find a point from a through b with distance d
PT point_along_line(PT a, PT b, double d) {
return a + (((b - a) / (b - a).norm()) * d);
}
// projection point c onto line through a and b assuming a != b
PT project_from_point_to_line(PT a, PT b, PT c) {
return a + (b - a) * dot(c - a, b - a) / (b - a).norm2();
}
// reflection point c onto line through a and b assuming a != b
PT reflection_from_point_to_line(PT a, PT b, PT c) {
PT p = project_from_point_to_line(a,b,c);
return point_along_line(c, p, 2.0 * dist(c, p));
}
// minimum distance from point c to line through a and b
double dist_from_point_to_line(PT a, PT b, PT c) {
return fabs(cross(b - a, c - a) / (b - a).norm());
}
// returns true if point p is on line segment ab
bool is_point_on_seg(PT a, PT b, PT p) {
if (fabs(cross(p - b, a - b)) < eps) {
if (p.x < min(a.x, b.x) || p.x > max(a.x, b.x)) return false;
if (p.y < min(a.y, b.y) || p.y > max(a.y, b.y)) return false;
return true;
}
return false;
}
// minimum distance point from point c to segment ab that lies on segment ab
PT project_from_point_to_seg(PT a, PT b, PT c) {
double r = dist2(a, b);
if (fabs(r) < eps) return a;
r = dot(c - a, b - a) / r;
if (r < 0) return a;
if (r > 1) return b;
return a + (b - a) * r;
}
// minimum distance from point c to segment ab
double dist_from_point_to_seg(PT a, PT b, PT c) {
return dist(c, project_from_point_to_seg(a, b, c));
}
// 0 if not parallel, 1 if parallel, 2 if collinear
bool is_parallel(PT a, PT b, PT c, PT d) {
double k = fabs(cross(b - a, d - c));
if (k < eps){
if (fabs(cross(a - b, a - c)) < eps && fabs(cross(c - d, c - a)) < eps) return 2;
else return 1;
}
else return 0;
}
// check if two lines are same
bool are_lines_same(PT a, PT b, PT c, PT d) {
if (fabs(cross(a - c, c - d)) < eps && fabs(cross(b - c, c - d)) < eps) return true;
return false;
}
// bisector vector of <abc
PT angle_bisector(PT &a, PT &b, PT &c){
PT p = a - b, q = c - b;
return p + q * sqrt(dot(p, p) / dot(q, q));
}
// 1 if point is ccw to the line, 2 if point is cw to the line, 3 if point is on the line
int point_line_relation(PT a, PT b, PT p) {
int c = sign(cross(p - a, b - a));
if (c < 0) return 1;
if (c > 0) return 2;
return 3;
}
// intersection point between ab and cd assuming unique intersection exists
bool line_line_intersection(PT a, PT b, PT c, PT d, PT &ans) {
double a1 = a.y - b.y, b1 = b.x - a.x, c1 = cross(a, b);
double a2 = c.y - d.y, b2 = d.x - c.x, c2 = cross(c, d);
double det = a1 * b2 - a2 * b1;
if (det == 0) return 0;
ans = PT((b1 * c2 - b2 * c1) / det, (c1 * a2 - a1 * c2) / det);
return 1;
}
// intersection point between segment ab and segment cd assuming unique intersection exists
bool seg_seg_intersection(PT a, PT b, PT c, PT d, PT &ans) {
double oa = cross2(c, d, a), ob = cross2(c, d, b);
double oc = cross2(a, b, c), od = cross2(a, b, d);
if (oa * ob < 0 && oc * od < 0){
ans = (a * ob - b * oa) / (ob - oa);
return 1;
}
else return 0;
}
// intersection point between segment ab and segment cd assuming unique intersection may not exists
// se.size()==0 means no intersection
// se.size()==1 means one intersection
// se.size()==2 means range intersection
set<PT> seg_seg_intersection_inside(PT a, PT b, PT c, PT d) {
PT ans;
if (seg_seg_intersection(a, b, c, d, ans)) return {ans};
set<PT> se;
if (is_point_on_seg(c, d, a)) se.insert(a);
if (is_point_on_seg(c, d, b)) se.insert(b);
if (is_point_on_seg(a, b, c)) se.insert(c);
if (is_point_on_seg(a, b, d)) se.insert(d);
return se;
}
// intersection between segment ab and line cd
// 0 if do not intersect, 1 if proper intersect, 2 if segment intersect
int seg_line_relation(PT a, PT b, PT c, PT d) {
double p = cross2(c, d, a);
double q = cross2(c, d, b);
if (sign(p) == 0 && sign(q) == 0) return 2;
else if (p * q < 0) return 1;
else return 0;
}
// intersection between segament ab and line cd assuming unique intersection exists
bool seg_line_intersection(PT a, PT b, PT c, PT d, PT &ans) {
bool k = seg_line_relation(a, b, c, d);
assert(k != 2);
if (k) line_line_intersection(a, b, c, d, ans);
return k;
}
// minimum distance from segment ab to segment cd
double dist_from_seg_to_seg(PT a, PT b, PT c, PT d) {
PT dummy;
if (seg_seg_intersection(a, b, c, d, dummy)) return 0.0;
else return min({dist_from_point_to_seg(a, b, c), dist_from_point_to_seg(a, b, d),
dist_from_point_to_seg(c, d, a), dist_from_point_to_seg(c, d, b)});
}
// minimum distance from point c to ray (starting point a and direction vector b)
double dist_from_point_to_ray(PT a, PT b, PT c) {
b = a + b;
double r = dot(c - a, b - a);
if (r < 0.0) return dist(c, a);
return dist_from_point_to_line(a, b, c);
}
// starting point as and direction vector ad
bool ray_ray_intersection(PT as, PT ad, PT bs, PT bd) {
double dx = bs.x - as.x, dy = bs.y - as.y;
double det = bd.x * ad.y - bd.y * ad.x;
if (fabs(det) < eps) return 0;
double u = (dy * bd.x - dx * bd.y) / det;
double v = (dy * ad.x - dx * ad.y) / det;
if (sign(u) >= 0 && sign(v) >= 0) return 1;
else return 0;
}
double ray_ray_distance(PT as, PT ad, PT bs, PT bd) {
if (ray_ray_intersection(as, ad, bs, bd)) return 0.0;
double ans = dist_from_point_to_ray(as, ad, bs);
ans = min(ans, dist_from_point_to_ray(bs, bd, as));
return ans;
}
60 changes: 60 additions & 0 deletions geometry/point.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
const double inf = 1e100;
const double eps = 1e-9;
const double PI = acos((double)-1.0);
int sign(double x) { return (x > eps) - (x < -eps); }
struct PT {
double x, y;
PT() { x = 0, y = 0; }
PT(double x, double y) : x(x), y(y) {}
PT(const PT &p) : x(p.x), y(p.y) {}
PT operator + (const PT &a) const { return PT(x + a.x, y + a.y); }
PT operator - (const PT &a) const { return PT(x - a.x, y - a.y); }
PT operator * (const double a) const { return PT(x * a, y * a); }
friend PT operator * (const double &a, const PT &b) { return PT(a * b.x, a * b.y); }
PT operator / (const double a) const { return PT(x / a, y / a); }
bool operator == (PT a) const { return sign(a.x - x) == 0 && sign(a.y - y) == 0; }
bool operator != (PT a) const { return !(*this == a); }
bool operator < (PT a) const { return sign(a.x - x) == 0 ? y < a.y : x < a.x; }
bool operator > (PT a) const { return sign(a.x - x) == 0 ? y > a.y : x > a.x; }
double norm() { return sqrt(x * x + y * y); }
double norm2() { return x * x + y * y; }
PT perp() { return PT(-y, x); }
double arg() { return atan2(y, x); }
PT truncate(double r) { // returns a vector with norm r and having same direction
double k = norm();
if (!sign(k)) return *this;
r /= k;
return PT(x * r, y * r);
}
};
inline double dot(PT a, PT b) { return a.x * b.x + a.y * b.y; }
inline double dist2(PT a, PT b) { return dot(a - b, a - b); }
inline double dist(PT a, PT b) { return sqrt(dot(a - b, a - b)); }
inline double cross(PT a, PT b) { return a.x * b.y - a.y * b.x; }
inline double cross2(PT a, PT b, PT c) { return cross(b - a, c - a); }
inline int orientation(PT a, PT b, PT c) { return sign(cross(b - a, c - a)); }
PT perp(PT a) { return PT(-a.y, a.x); }
PT rotateccw90(PT a) { return PT(-a.y, a.x); }
PT rotatecw90(PT a) { return PT(a.y, -a.x); }
PT rotateccw(PT a, double t) { return PT(a.x * cos(t) - a.y * sin(t), a.x * sin(t) + a.y * cos(t)); }
PT rotatecw(PT a, double t) { return PT(a.x * cos(t) + a.y * sin(t), -a.x * sin(t) + a.y * cos(t)); }
double SQ(double x) { return x * x; }
double rad_to_deg(double r) { return (r * 180.0 / PI); }
double deg_to_rad(double d) { return (d * PI / 180.0); }
double get_angle(PT a, PT b) {
double costheta = dot(a, b) / a.norm() / b.norm();
return acos(max((double)-1.0, min((double)1.0, costheta)));
}
bool is_point_in_angle(PT b, PT a, PT c, PT p) { // does point p lie in angle <bac
assert(orientation(a, b, c) != 0);
if (orientation(a, c, b) < 0) swap(b, c);
return orientation(a, c, p) >= 0 && orientation(a, b, p) <= 0;
}
bool half(PT p) {
return p.y > 0.0 || (p.y == 0.0 && p.x < 0.0);
}
void polar_sort(vector<PT> &v) { // sort points in counterclockwise
sort(v.begin(), v.end(), [](PT a,PT b) {
return make_tuple(half(a), 0.0, a.norm2()) < make_tuple(half(b), cross(a, b), b.norm2());
});
}
Loading

0 comments on commit 14ea361

Please sign in to comment.