furthest_pair.hpp

SECTIONGeometry INCLUDEnoya/furthest_pair.hpp

求平面点集中欧氏距离最远的一对点；通常在凸包上用旋转卡壳完成。

Complexity: Time: O(n log n). Space: O(n).

AC 记录：furthest_pair。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @complexity Time: O(n log n).
/// Space: O(n).

#include "noya/geometry_base.hpp"

#include <algorithm>
#include <cassert>
#include <optional>
#include <utility>
#include <vector>

namespace noya {

/// @brief Return indices of a furthest pair in O(n log n).
template <class T>
std::pair<int, int> furthest_pair(const std::vector<point<T>> &pt) {
  assert(pt.size() >= 2);
  std::vector<point<T>> h = convex_hull(pt);
  if (h.size() == 1) {
    return {0, 1};
  }

  std::vector<std::pair<point<T>, int>> ids;
  ids.reserve(pt.size());
  for (int i = 0; i < int(pt.size()); i++) {
    ids.emplace_back(pt[i], i);
  }
  std::sort(ids.begin(), ids.end());
  std::vector<int> idx;
  idx.reserve(h.size());
  for (const point<T> &val : h) {
    auto it = std::lower_bound(
        ids.begin(), ids.end(), val,
        [](const auto &ent, const point<T> &can) { return ent.first < can; });
    idx.push_back(it->second);
  }

  using distance_type = decltype(norm2(point<T>{}));
  std::optional<distance_type> bst;
  std::pair<int, int> ans = {-1, -1};
  auto upd = [&](int a, int b) {
    distance_type can = norm2(h[a] - h[b]);
    std::pair<int, int> id1 = std::minmax(idx[a], idx[b]);
    if (!bst || can > *bst || (can == *bst && id1 < ans)) {
      bst = can;
      ans = id1;
    }
  };

  const int N = int(h.size());
  if (N == 2) {
    upd(0, 1);
    return ans;
  }
  int opp = 1;
  for (int i = 0; i < N; i++) {
    int ni = (i + 1) % N;
    while (true) {
      int nj = (opp + 1) % N;
      auto s = cross(h[ni] - h[i], h[opp] - h[i]);
      auto na = cross(h[ni] - h[i], h[nj] - h[i]);
      if (na <= s) {
        break;
      }
      opp = nj;
    }
    upd(i, opp);
    upd(ni, opp);
  }
  return ans;
}

} // namespace noya

#ifndef NOYA_FURTHEST_PAIR_HPP
#define NOYA_FURTHEST_PAIR_HPP 1

/// @complexity Time: O(n log n).
/// Space: O(n).

#include "noya/geometry_base.hpp"

#include <algorithm>
#include <cassert>
#include <optional>
#include <utility>
#include <vector>

namespace noya {

/// @brief Return indices of a furthest pair in O(n log n).
template <class T>
std::pair<int, int> furthest_pair(const std::vector<point<T>> &pt) {
  assert(pt.size() >= 2);
  std::vector<point<T>> h = convex_hull(pt);
  if (h.size() == 1) {
    return {0, 1};
  }

  std::vector<std::pair<point<T>, int>> ids;
  ids.reserve(pt.size());
  for (int i = 0; i < int(pt.size()); i++) {
    ids.emplace_back(pt[i], i);
  }
  std::sort(ids.begin(), ids.end());
  std::vector<int> idx;
  idx.reserve(h.size());
  for (const point<T> &val : h) {
    auto it = std::lower_bound(
        ids.begin(), ids.end(), val,
        [](const auto &ent, const point<T> &can) { return ent.first < can; });
    idx.push_back(it->second);
  }

  using distance_type = decltype(norm2(point<T>{}));
  std::optional<distance_type> bst;
  std::pair<int, int> ans = {-1, -1};
  auto upd = [&](int a, int b) {
    distance_type can = norm2(h[a] - h[b]);
    std::pair<int, int> id1 = std::minmax(idx[a], idx[b]);
    if (!bst || can > *bst || (can == *bst && id1 < ans)) {
      bst = can;
      ans = id1;
    }
  };

  const int N = int(h.size());
  if (N == 2) {
    upd(0, 1);
    return ans;
  }
  int opp = 1;
  for (int i = 0; i < N; i++) {
    int ni = (i + 1) % N;
    while (true) {
      int nj = (opp + 1) % N;
      auto s = cross(h[ni] - h[i], h[opp] - h[i]);
      auto na = cross(h[ni] - h[i], h[nj] - h[i]);
      if (na <= s) {
        break;
      }
      opp = nj;
    }
    upd(i, opp);
    upd(ni, opp);
  }
  return ans;
}

} // namespace noya

#endif // NOYA_FURTHEST_PAIR_HPP

#include <algorithm>
#include <cassert>
#include <cmath>
#include <optional>
#include <utility>
#include <vector>

/// @complexity Time: O(n log n).
/// Space: O(n).

/// @complexity Time: O(1) per primitive; O(n log n) for convex hull.
/// Space: O(1) per primitive and O(n) for hull construction.

namespace noya {

/// @brief Two-dimensional point with vector arithmetic and lexicographic order.
template <class T> struct point {
  T x{};
  T y{};

  point() = default;
  point(T x_, T y_) : x(x_), y(y_) {}

  point &operator+=(const point &rhs) {
    x += rhs.x;
    y += rhs.y;
    return *this;
  }
  point &operator-=(const point &rhs) {
    x -= rhs.x;
    y -= rhs.y;
    return *this;
  }
  point &operator*=(const T &scl) {
    x *= scl;
    y *= scl;
    return *this;
  }
  point &operator/=(const T &scl) {
    x /= scl;
    y /= scl;
    return *this;
  }

  friend point operator+(point l, const point &r) { return l += r; }
  friend point operator-(point l, const point &r) { return l -= r; }
  friend point operator*(point val, const T &scl) { return val *= scl; }
  friend point operator*(const T &scl, point val) { return val *= scl; }
  friend point operator/(point val, const T &scl) { return val /= scl; }
  friend bool operator==(const point &, const point &) = default;
  friend bool operator<(const point &l, const point &r) {
    return l.x < r.x || (l.x == r.x && l.y < r.y);
  }
};

/// @brief Circle represented by a center and a nonnegative radius.
template <class Real> struct circle {
  point<Real> o;
  Real r1{};
};

/// @brief Return the dot product of two vectors.
template <class T> T dot(const point<T> &a, const point<T> &b) {
  return a.x * b.x + a.y * b.y;
}

/// @brief Return the signed cross product of two vectors.
template <class T> T cross(const point<T> &a, const point<T> &b) {
  return a.x * b.y - a.y * b.x;
}

/// @brief Return cross(a - o1, b - o1).
template <class T>
T cross(const point<T> &o1, const point<T> &a, const point<T> &b) {
  return cross(a - o1, b - o1);
}

/// @brief Return the squared Euclidean norm.
template <class T> T norm2(const point<T> &val) { return dot(val, val); }

/// @brief Compare a value with zero using an optional absolute tolerance.
template <class T> int sign(const T &val, const T &eps = T{}) {
  return (val > eps) - (val < -eps);
}

/// @brief Return -1, 0, or 1 for a clockwise, collinear, or counter-clockwise
/// turn.
template <class T>
int orientation(const point<T> &a, const point<T> &b, const point<T> &c,
                const T &eps = T{}) {
  return sign(cross(a, b, c), eps);
}

/// @brief Test whether p lies on the closed segment [a, b].
template <class T>
bool on_segment(const point<T> &p, const point<T> &a, const point<T> &b,
                const T &eps = T{}) {
  if (orientation(a, b, p, eps) != 0) {
    return false;
  }
  return std::min(a.x, b.x) - eps <= p.x && p.x <= std::max(a.x, b.x) + eps &&
         std::min(a.y, b.y) - eps <= p.y && p.y <= std::max(a.y, b.y) + eps;
}

/// @brief Test whether the closed segments [a, b] and [c, d] intersect.
template <class T>
bool segments_intersect(const point<T> &a, const point<T> &b, const point<T> &c,
                        const point<T> &d, const T &eps = T{}) {
  int s1 = orientation(a, b, c, eps);
  int s2 = orientation(a, b, d, eps);
  int s3 = orientation(c, d, a, eps);
  int s4 = orientation(c, d, b, eps);
  if (s1 == 0 && on_segment(c, a, b, eps)) {
    return true;
  }
  if (s2 == 0 && on_segment(d, a, b, eps)) {
    return true;
  }
  if (s3 == 0 && on_segment(a, c, d, eps)) {
    return true;
  }
  if (s4 == 0 && on_segment(b, c, d, eps)) {
    return true;
  }
  return s1 * s2 < 0 && s3 * s4 < 0;
}

/// @brief Intersect the infinite lines through (a, b) and (c, d), returning
/// nullopt when they are parallel or coincident.
template <class T>
std::optional<point<long double>>
line_intersection(const point<T> &a, const point<T> &b, const point<T> &c,
                  const point<T> &d, long double eps = 0) {
  point<long double> lhs{static_cast<long double>(a.x),
                         static_cast<long double>(a.y)};
  point<long double> b1{static_cast<long double>(b.x),
                        static_cast<long double>(b.y)};
  point<long double> z{static_cast<long double>(c.x),
                       static_cast<long double>(c.y)};
  point<long double> d1{static_cast<long double>(d.x),
                        static_cast<long double>(d.y)};
  point<long double> da = b1 - lhs;
  point<long double> db = d1 - z;
  long double den = cross(da, db);
  if (std::abs(den) <= eps) {
    return std::nullopt;
  }
  long double rat = cross(z - lhs, db) / den;
  return lhs + da * rat;
}

/// @brief Return the convex hull in counter-clockwise order without repetition.
template <class T>
std::vector<point<T>> convex_hull(std::vector<point<T>> pt, bool kc = false) {
  std::sort(pt.begin(), pt.end());
  pt.erase(std::unique(pt.begin(), pt.end()), pt.end());
  if (pt.size() <= 1) {
    return pt;
  }
  bool col = true;
  for (int i = 2; i < int(pt.size()); i++) {
    col &= orientation(pt[0], pt[1], pt[i]) == 0;
  }
  if (kc && col) {
    return pt;
  }
  std::vector<point<T>> lo, hi;
  for (const point<T> &p : pt) {
    while (lo.size() >= 2) {
      int rot = orientation(lo[lo.size() - 2], lo.back(), p);
      if (rot > 0 || (kc && rot == 0)) {
        break;
      }
      lo.pop_back();
    }
    lo.push_back(p);
  }
  for (auto it = pt.rbegin(); it != pt.rend(); ++it) {
    while (hi.size() >= 2) {
      int rot = orientation(hi[hi.size() - 2], hi.back(), *it);
      if (rot > 0 || (kc && rot == 0)) {
        break;
      }
      hi.pop_back();
    }
    hi.push_back(*it);
  }
  lo.pop_back();
  hi.pop_back();
  lo.insert(lo.end(), hi.begin(), hi.end());
  return lo;
}

/// @brief Return twice the signed area of a polygon.
template <class T> T polygon_area2(const std::vector<point<T>> &pg) {
  T res{};
  for (int i = 0; i < int(pg.size()); i++) {
    res += cross(pg[i], pg[(i + 1) % pg.size()]);
  }
  return res;
}

/// @brief Classify a point relative to a polygon: -1 outside, 0 boundary, 1
/// inside.
template <class T>
int point_in_polygon(const point<T> &p, const std::vector<point<T>> &pg) {
  bool in = false;
  for (int i = 0; i < int(pg.size()); i++) {
    point<T> a = pg[i];
    point<T> b = pg[(i + 1) % pg.size()];
    if (on_segment(p, a, b)) {
      return 0;
    }
    if (a.y <= p.y && p.y < b.y && orientation(a, b, p) > 0) {
      in = !in;
    }
    if (b.y <= p.y && p.y < a.y && orientation(a, b, p) < 0) {
      in = !in;
    }
  }
  return in ? 1 : -1;
}

/// @brief Return the Euclidean distance between two points.
template <class T> long double distance(const point<T> &a, const point<T> &b) {
  return std::hypot(
      static_cast<long double>(a.x) - static_cast<long double>(b.x),
      static_cast<long double>(a.y) - static_cast<long double>(b.y));
}

} // namespace noya

namespace noya {

/// @brief Return indices of a furthest pair in O(n log n).
template <class T>
std::pair<int, int> furthest_pair(const std::vector<point<T>> &pt) {
  assert(pt.size() >= 2);
  std::vector<point<T>> h = convex_hull(pt);
  if (h.size() == 1) {
    return {0, 1};
  }

  std::vector<std::pair<point<T>, int>> ids;
  ids.reserve(pt.size());
  for (int i = 0; i < int(pt.size()); i++) {
    ids.emplace_back(pt[i], i);
  }
  std::sort(ids.begin(), ids.end());
  std::vector<int> idx;
  idx.reserve(h.size());
  for (const point<T> &val : h) {
    auto it = std::lower_bound(
        ids.begin(), ids.end(), val,
        [](const auto &ent, const point<T> &can) { return ent.first < can; });
    idx.push_back(it->second);
  }

  using distance_type = decltype(norm2(point<T>{}));
  std::optional<distance_type> bst;
  std::pair<int, int> ans = {-1, -1};
  auto upd = [&](int a, int b) {
    distance_type can = norm2(h[a] - h[b]);
    std::pair<int, int> id1 = std::minmax(idx[a], idx[b]);
    if (!bst || can > *bst || (can == *bst && id1 < ans)) {
      bst = can;
      ans = id1;
    }
  };

  const int N = int(h.size());
  if (N == 2) {
    upd(0, 1);
    return ans;
  }
  int opp = 1;
  for (int i = 0; i < N; i++) {
    int ni = (i + 1) % N;
    while (true) {
      int nj = (opp + 1) % N;
      auto s = cross(h[ni] - h[i], h[opp] - h[i]);
      auto na = cross(h[ni] - h[i], h[nj] - h[i]);
      if (na <= s) {
        break;
      }
      opp = nj;
    }
    upd(i, opp);
    upd(ni, opp);
  }
  return ans;
}

} // namespace noya