【计算几何】全代码

const double eps = 1e-7;
const double pi = acos(-1);
int dcmp(double x)
{ // 判负
  return (fabs(x) <= eps) ? 0 : (x < 0 ? -1 : 1);
}
struct Point
{ // 点
  double x, y;
  Point(double X = 0, double Y = 0) { x = X, y = Y; }
};
struct Vector
{ // 向量
  double x, y;
  Vector(double X = 0, double Y = 0) { x = X, y = Y; }
};
inline Vector operator+(Vector x, Vector y)
{ // 向量+向量=向量
  return Vector(x.x + y.x, x.y + y.y);
}
inline Point operator+(Point x, Vector y)
{ // 点+向量=点
  return Point(x.x + y.x, x.y + y.y);
}
inline Point operator+(Vector y, Point x)
{ // 向量+点=点
  return Point(x.x + y.x, x.y + y.y);
}
inline Vector operator-(Vector x, Vector y)
{ // 向量-向量=向量
  return Vector(x.x - y.x, x.y - y.y);
}
inline Point operator-(Point x, Vector y)
{ // 点-向量=点
  return Point(x.x - y.x, x.y - y.y);
}
inline Vector operator-(Point x, Point y)
{ // 点-点=向量
  return Vector(x.x - y.x, x.y - y.y);
}
inline Vector operator*(Vector x, double y)
{ // 向量*数=向量
  return Vector(x.x * y, x.y * y);
}
inline Vector operator/(Vector x, double y)
{ // 向量/数=向量
  return Vector(x.x / y, x.y / y);
}
inline double polar(Vector x)
{ // 向量极角
  return atan2(x.y, x.x);
}
inline double len(Vector x)
{ // 向量模长
  return sqrt(x.x * x.x + x.y * x.y);
}
inline double dot(Vector x, Vector y)
{ // 向量点积
  return x.x * y.x + x.y * y.y;
}
inline double cross(Vector x, Vector y)
{ // 向量叉积
  return x.x * y.y - x.y * y.x;
}
inline double operator*(Vector x, Vector y)
{ // 向量叉积
  return cross(x, y);
}
inline double theta(Vector x, Vector y)
{ // 向量夹角
  return acos(dot(x, y) / len(x) / len(y));
}
inline Vector rotate(Vector A, double rad)
{ // 向量旋转,rad 是顺时针弧度制
  return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));
}
struct Line
{
  Point x;
  Vector y;
  Line(Point X, Vector Y) { x = X, y = Y; }
  Line(Point X, Point Y) { x = X, y = Y - X; }
};
inline Point GetLineIntersection(Line x, Line y)
{ // 两直线交点
  Vector u = x.x - y.x;
  double t = cross(y.y, u) / cross(x.y, y.y);
  return x.x + x.y * t;
}
inline double DistanceToLine(Point P, Line x)
{ // 点到直线的距离
  Vector v1 = x.y, v2 = P - x.x;
  return fabs(cross(v1, v2)) / len(v1);
}
inline Point GetLineProjection(Point P, Line x)
{ // 点到直线的投影
  return x.x + x.y * (dot(x.y, P - x.x) / dot(x.y, x.y));
}
double PolygonSignedArea(vector<Point> poly)
{ // 多边形有向面积
  double area = 0;
  for (int i = 1; i < (poly.size() - 1); i++)
  {
    area += (poly[i] - poly[0]) * (poly[i + 1] - poly[0]);
  }
  return area / 2;
}
double PolygonArea(vector<Point> poly)
{ // 多边形无向面积
  double area = 0;
  for (int i = 1; i < (poly.size() - 1); i++)
  {
    area += fabs((poly[i] - poly[0]) * (poly[i + 1] - poly[0]));
  }
  return area / 2;
}
inline int Andrew(Point *p, int n, Point *Poly) { // 求凸包
  int used[N] = {0};                              // 标记每个数字是否使用过
  sort(p + 1, p + n + 1, [](Point a, Point b) {   // 排序节点
    return a.x < b.x || (a.x == b.x && a.y < b.y);
  });
  int cnt = 3;
  Poly[1] = p[1], Poly[2] = p[2];                 // 初始时加入1号和2号节点
  for (int i = 3; i <= n; ++i) {
    if (cnt > 1 && Cross(Poly[cnt] - Poly[cnt - 1], Poly[i] - Poly[cnt] <= 0)) { // 判断是否在左边
      -- cnt;                                                                    // 删除此节点
      used[Poly[cnt + 1]] = 0;
    }
    used[i] = 1;                                 // 标记节点
    Poly[++ cnt] = p[i];                         // 加入凸包
  }
  int k = cnt;                                   // 上壳的起点
  for (int i = n; i >= 1; --i) {
    if (used[i]) {                               // 已经加入凸包
      continue;
    }
    while (cnt > k && Cross(Poly[cnt] - Poly[cnt - 1], Poly[i] - Poly[cnt]) <= 0)) {
      -- cnt;
      used[Poly[cnt + 1]] = 0;
    }
    used[i] = 1;
    Poly[++ cnt] = p[i];
  }
  Poly[++ cnt] = 1;                              // 便于计算周长
  return cnt;
}

posted @ 2023-01-12 09:11  SenGYi  阅读(32)  评论(0)    收藏  举报