计算几何

hj23308 · 2024-04-17 14:59:53 · 算法·理论

basic

template <typename D> int Sign(D x) {
  constexpr D EPS = std::is_integral<D>::value ? 0 : 1E-10;
  /*
    if (x > +EPS) return +1;
    if (x < -EPS) return -1;
    return 0;
  */
  return (x > +EPS) - (x < -EPS);
}

template<typename D>
struct Vector {
  D x{}, y{};

  Vector() = default;
  Vector(D x, D y) : x(x), y(y) {}

  Vector operator-() const { return {-x, -y}; }

  Vector operator*(D scalar) const { return {x * scalar, y * scalar}; }
  Vector operator/(D scalar) const { return {x / scalar, y / scalar}; }

  Vector operator+(const Vector& rhs) const { return {x + rhs.x, y + rhs.y}; }
  Vector operator-(const Vector& rhs) const { return {x - rhs.x, y - rhs.y}; }

  friend D dot(Vector a, Vector b) { return a.x * b.x + a.y * b.y; }
  friend D crs(Vector a, Vector b) { return a.x * b.y - a.y * b.x; }

  D sqd() const { return x * x + y * y; }
  D len() const { return std::sqrt(sqd()); }
  D ang() const { return std::atan2(y, x); }

  Vector norm() const { return (*this) / len(); }

  friend bool par(Vector a, Vector b) { return Sign(crs(a, b)) == 0; }

  friend std::ostream& operator<<(std::ostream& os, const Vector& v) {
    return os << '(' << v.x << ", " << v.y << ')';
  }
};

template<typename D>
struct Segment {
  D angle;
  Vector<D> s, t;

  Segment() = default;
  Segment(Vector<D> s, Vector<D> t) : s(s), t(t) {}

  D ang() { return angle = (t - s).ang(); }

  int side(Vector<D> p) { return Sign(crs(p - s, t - s)); }

  friend bool par(Segment a, Segment b) { return par(a.t - a.s, b.t - b.s); }

  friend Vector<D> its(Segment a, Segment b) {
    D k = crs(b.t - b.s, a.s - b.s) / crs(a.t - a.s, b.t - b.s);
    return a.s + (a.t - a.s) * k;
  }

  friend std::ostream& operator<<(std::ostream& os, const Segment& p) {
    return os << '{' << p.s << ", " << p.t << '}';
  }
};

polygonal area

template<typename D>
D Area(const std::vector<Vector<D>>& p) {
  D s = 0;
  int n = (int) p.size();
  for (int i = 0; i < n; i++)
    s += crs(p[i], p[(i + 1) % n]);
  return s / 2;
}

build convex hull

template<typename D>
std::vector<Vector<D>> ConvexHull(std::vector<Vector<D>> points) {
  std::sort(points.begin(), points.end(), [&](auto a, auto b) {
    return std::tie(a.x, a.y) < std::tie(b.x, b.y);
  });

  std::vector<Vector<D>> result;
  for (int step = 0; step < 2; step++) {
    int top = 0;
    std::vector<Vector<D>> sta(points.size());
    for (const auto& p : points) {
      while (top > 1 && Sign(crs(sta[top - 2] - sta[top - 1],
        sta[top - 2] - p)) < 1) top--;
      sta[top++] = p;
    }
    for (int i = 0; i < top - 1; i++) result.emplace_back(sta[i]);
    std::reverse(points.begin(), points.end());
  }
  return result;
}

polar angle sorting

template<typename D>
inline void PASort(std::vector<Vector<D>>& points) {
  auto quad = [&](Vector<D> p) {
    return Sign(p.y) == 1 || (!Sign(p.y) && Sign(p.x) < 0);
  };
  std::sort(points.begin(), points.end(), [&](auto a, auto b) {
    return quad(a) == quad(b) ? Sign(crs(a, b)) > 0 : quad(b);
  });
}

half plane intersection

template<typename D>
std::vector<Vector<D>> HPIts(std::vector<Segment<D>> segs, bool ext = true) {
  if (ext) {
    constexpr D INF = 1E18;
    Vector<D> cur(INF, INF);
    for (int step = 0; step < 4; step++) {
      Vector<D> nxt(-cur.y, cur.x);
      segs.emplace_back(cur, nxt);
      cur = nxt;
    }
  }

  for (auto& s : segs) s.ang();
  std::sort(segs.begin(), segs.end(), [&](auto a, auto b) {
    return Sign(a.angle - b.angle) ? a.angle < b.angle : a.side(b.s) == 1;
  });

  std::deque<Vector<D>> points;
  std::deque<Segment<D>> queue;
  std::vector<Vector<D>> result;

  for (int i = 0; i < (int) segs.size(); i++) {
    if (i && !Sign(segs[i].angle - segs[i - 1].angle))
      continue;
    while (queue.size() > 1 && segs[i].side(points.back()) == 1)
      queue.pop_back(), points.pop_back();
    while (queue.size() > 1 && segs[i].side(points.front()) == 1)
      queue.pop_front(), points.pop_front();
    if (i) {
      if (par(queue.back(), segs[i])) return result;
      points.emplace_back(its(queue.back(), segs[i]));
    }
    queue.emplace_back(segs[i]);
  }
  while (queue.size() > 1 && queue.front().side(points.back()) == 1)
    queue.pop_back(), points.pop_back();

  if (queue.size() > 2) {
    for (auto v : points)
      result.emplace_back(v);
    result.emplace_back(its(queue.front(), queue.back()));
  }
  return result;
}