geometry/segment.hpp

• View this file on GitHub
• Last update: 2024-05-03 21:06:15+09:00
• Include: #include "geometry/segment.hpp"

Depends on

• geometry/geometry-base.hpp
• geometry/line.hpp
• geometry/polygon.hpp

Verified with

• verify/verify-aoj-cgl/aoj-cgl-2-b.test.cpp
• verify/verify-aoj-cgl/aoj-cgl-2-d.test.cpp

Code

#pragma once

#include "geometry-base.hpp"
//
#include "line.hpp"

struct Segment : Line {
  Segment() = default;
  using Line::Line;
};

using Segments = vector<Segment>;

bool is_intersect_sp(const Segment &s, const Point &p) {
  return ccw(s.a, s.b, p) == 0;
}
bool is_intersect_ss(const Segment &s, const Segment &t) {
  return ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) <= 0 &&
         ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) <= 0;
}

Real distance_sp(const Segment &s, const Point &p) {
  Point r = projection(s, p);
  if (is_intersect_sp(s, r)) return abs(r - p);
  return min(abs(s.a - p), abs(s.b - p));
}
Real distance_ss(const Segment &a, const Segment &b) {
  if (is_intersect_ss(a, b)) return 0;
  return min({distance_sp(a, b.a), distance_sp(a, b.b), distance_sp(b, a.a),
              distance_sp(b, a.b)});
}

#line 2 "geometry/segment.hpp"

#line 2 "geometry/geometry-base.hpp"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <complex>
#include <iostream>
#include <vector>
using namespace std;

using Real = long double;
constexpr Real EPS = 1e-10;
constexpr Real pi = 3.141592653589793238462643383279L;
bool equals(Real a, Real b) { return fabs(b - a) < EPS; }
int sign(Real a) { return equals(a, 0) ? 0 : a > 0 ? 1 : -1; }

template <typename R>
struct PointBase {
  using P = PointBase;
  R x, y;
  PointBase() : x(0), y(0) {}
  PointBase(R _x, R _y) : x(_x), y(_y) {}
  template <typename T, typename U>
  PointBase(const pair<T, U>& p) : x(p.first), y(p.second) {}

  P operator+(const P& r) const { return {x + r.x, y + r.y}; }
  P operator-(const P& r) const { return {x - r.x, y - r.y}; }
  P operator*(R r) const { return {x * r, y * r}; }
  P operator/(R r) const { return {x / r, y / r}; }

  P& operator+=(const P& r) { return (*this) = (*this) + r; }
  P& operator-=(const P& r) { return (*this) = (*this) - r; }
  P& operator*=(R r) { return (*this) = (*this) * r; }
  P& operator/=(R r) { return (*this) = (*this) / r; }

  bool operator<(const P& r) const { return x != r.x ? x < r.x : y < r.y; }
  bool operator==(const P& r) const { return x == r.x and y == r.y; }
  bool operator!=(const P& r) const { return !((*this) == r); }

  P rotate(R rad) const {
    return {x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad)};
  }

  R real() const { return x; }
  R imag() const { return y; }
  friend R real(const P& p) { return p.x; }
  friend R imag(const P& p) { return p.y; }
  friend R dot(const P& l, const P& r) { return l.x * r.x + l.y * r.y; }
  friend R cross(const P& l, const P& r) { return l.x * r.y - l.y * r.x; }
  friend R abs(const P& p) { return sqrt(p.x * p.x + p.y * p.y); }
  friend R norm(const P& p) { return p.x * p.x + p.y * p.y; }
  friend R arg(const P& p) { return atan2(p.y, p.x); }

  friend istream& operator>>(istream& is, P& p) {
    R a, b;
    is >> a >> b;
    p = P{a, b};
    return is;
  }
  friend ostream& operator<<(ostream& os, const P& p) {
    return os << p.x << " " << p.y;
  }
};
using Point = PointBase<Real>;
using Points = vector<Point>;

// ccw, 点の進行方向
int ccw(const Point& a, const Point& b, const Point& c) {
  Point x = b - a, y = c - a;
  if (cross(x, y) > EPS) return +1;                 // 反時計回り
  if (cross(x, y) < -EPS) return -1;                // 時計回り
  if (min(norm(x), norm(y)) < EPS * EPS) return 0;  // c=a または c=b
  if (dot(x, y) < EPS) return +2;                   // c-a-b の順で一直線
  if (norm(x) < norm(y)) return -2;                 // a-b-c の順で一直線
  return 0;                                         // a-c-b の順で一直線
}
#line 4 "geometry/segment.hpp"
//
#line 2 "geometry/line.hpp"

#line 2 "geometry/polygon.hpp"

#line 4 "geometry/polygon.hpp"

using Polygon = vector<Point>;

// 多角形の内部に点があるか？
// OUT : 0, ON : 1, IN : 2
int contains_polygon(const Polygon &Q, const Point &p) {
  bool in = false;
  for (int i = 0; i < (int)Q.size(); i++) {
    Point a = Q[i] - p, b = Q[(i + 1) % Q.size()] - p;
    if (imag(a) > imag(b)) swap(a, b);
    if (sign(imag(a)) <= 0 && 0 < sign(imag(b)) && sign(cross(a, b)) < 0)
      in = !in;
    if (equals(cross(a, b), 0) && sign(dot(a, b)) <= 0) return 1;
  }
  return in ? 2 : 0;
}

// 多角形の面積
Real area(const Polygon &p) {
  Real A = 0;
  for (int i = 0; i < (int)p.size(); ++i) {
    A += cross(p[i], p[(i + 1) % p.size()]);
  }
  return A * 0.5;
}

// 頂点集合から凸包を生成
// boundary : 周上の点も列挙する場合 true
template <bool boundary = false>
Polygon convex_hull(vector<Point> ps) {
  sort(begin(ps), end(ps));
  ps.erase(unique(begin(ps), end(ps)), end(ps));
  int n = ps.size(), k = 0;
  if (n <= 2) return ps;
  vector<Point> ch(2 * n);
  // 反時計周り
  Real th = boundary ? -EPS : +EPS;
  for (int i = 0; i < n; ch[k++] = ps[i++]) {
    while (k >= 2 && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k;
  }
  for (int i = n - 2, t = k + 1; i >= 0; ch[k++] = ps[i--]) {
    while (k >= t && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k;
  }
  ch.resize(k - 1);
  return ch;
}

// 凸包の内部に点があるか？
// OUT : 0, ON : 1, IN : 2
int contains_convex(const Polygon &C, const Point &p) {
  int N = C.size();
  auto b1 = cross(C[1] - C[0], p - C[0]);
  auto b2 = cross(C[N - 1] - C[0], p - C[0]);
  if (b1 < -EPS or b2 > EPS) return 0;
  int L = 1, R = N - 1;
  while (L + 1 < R) {
    int M = (L + R) / 2;
    (cross(p - C[0], C[M] - C[0]) >= 0 ? R : L) = M;
  }
  auto v = cross(C[L] - p, C[R] - p);
  if (equals(v, 0)) {
    return 1;
  } else if (v > 0) {
    return equals(b1, 0) or equals(b2, 0) ? 1 : 2;
  } else {
    return 0;
  }
}

// 凸包が与えられるので最遠点対を返す
// 返り値：頂点番号のペア
pair<int, int> convex_polygon_diameter(const Polygon &p) {
  int N = (int)p.size();
  int is = 0, js = 0;
  for (int i = 1; i < N; i++) {
    if (imag(p[i]) > imag(p[is])) is = i;
    if (imag(p[i]) < imag(p[js])) js = i;
  }
  Real maxdis = norm(p[is] - p[js]);

  int maxi, maxj, i, j;
  i = maxi = is;
  j = maxj = js;
  do {
    if (cross(p[(i + 1) % N] - p[i], p[(j + 1) % N] - p[j]) >= 0) {
      j = (j + 1) % N;
    } else {
      i = (i + 1) % N;
    }
    if (norm(p[i] - p[j]) > maxdis) {
      maxdis = norm(p[i] - p[j]);
      maxi = i;
      maxj = j;
    }
  } while (i != is || j != js);
  return minmax(maxi, maxj);
}
#line 5 "geometry/line.hpp"

struct Line {
  Point a, b;

  Line() = default;
  Line(const Point &_a, const Point &_b) : a(_a), b(_b) {}
  // Ax+By=C
  Line(const Real &A, const Real &B, const Real &C) {
    if (equals(A, 0)) {
      assert(!equals(B, 0));
      a = Point(0, C / B);
      b = Point(1, C / B);
    } else if (equals(B, 0)) {
      a = Point(C / A, 0);
      b = Point(C / A, 1);
    } else if (equals(C, 0)) {
      a = Point(0, C / B);
      b = Point(1, (C - A) / B);
    } else {
      a = Point(0, C / B);
      b = Point(C / A, 0);
    }
  }
  friend ostream &operator<<(ostream &os, const Line &l) {
    return os << l.a << " to " << l.b;
  }
  friend istream &operator>>(istream &is, Line &l) { return is >> l.a >> l.b; }
};
using Lines = vector<Line>;

bool is_parallel(const Line &a, const Line &b) {
  return equals(cross(a.b - a.a, b.b - b.a), 0);
}
bool is_orthogonal(const Line &a, const Line &b) {
  return equals(dot(a.a - a.b, b.a - b.b), 0);
}
Point cross_point_ll(const Line &l, const Line &m) {
  Real A = cross(l.b - l.a, m.b - m.a);
  Real B = cross(l.b - l.a, l.b - m.a);
  if (equals(abs(A), 0) && equals(abs(B), 0)) return m.a;
  return m.a + (m.b - m.a) * B / A;
}
bool is_intersect_ll(const Line &l, const Line &m) {
  Real A = cross(l.b - l.a, m.b - m.a);
  Real B = cross(l.b - l.a, l.b - m.a);
  if (equals(abs(A), 0) && equals(abs(B), 0)) return true;
  return !is_parallel(l, m);
}

// 直線に頂点から垂線を下ろした時の交点
Point projection(const Line &l, const Point &p) {
  Real t = dot(p - l.a, l.a - l.b) / norm(l.a - l.b);
  return l.a + (l.a - l.b) * t;
}

// 凸包を直線で切った時の片方 (直線 a->b の進行方向左側) を返す
Polygon convex_polygon_cut(const Polygon &U, const Line &l) {
  Polygon ret;
  for (int i = 0; i < (int)U.size(); i++) {
    const Point &now = U[i];
    const Point &nxt = U[(i + 1) % U.size()];
    auto cf = cross(l.a - now, l.b - now);
    auto cs = cross(l.a - nxt, l.b - nxt);
    if (sign(cf) >= 0) {
      ret.emplace_back(now);
    }
    if (sign(cf) * sign(cs) < 0) {
      ret.emplace_back(cross_point_ll(Line(now, nxt), l));
    }
  }
  return ret;
}
#line 6 "geometry/segment.hpp"

struct Segment : Line {
  Segment() = default;
  using Line::Line;
};

using Segments = vector<Segment>;

bool is_intersect_sp(const Segment &s, const Point &p) {
  return ccw(s.a, s.b, p) == 0;
}
bool is_intersect_ss(const Segment &s, const Segment &t) {
  return ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) <= 0 &&
         ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) <= 0;
}

Real distance_sp(const Segment &s, const Point &p) {
  Point r = projection(s, p);
  if (is_intersect_sp(s, r)) return abs(r - p);
  return min(abs(s.a - p), abs(s.b - p));
}
Real distance_ss(const Segment &a, const Segment &b) {
  if (is_intersect_ss(a, b)) return 0;
  return min({distance_sp(a, b.a), distance_sp(a, b.b), distance_sp(b, a.a),
              distance_sp(b, a.b)});
}

Back to top page