这是一套基础的计算几何模板组件, 可以解决的问题有
0. 内积和外积计算
1. 平行和正交判定.
2. 找投影点
3. 找对称点
4. 找距离(点与点, 点与直线, 点与线段, 线段与线段)
5. 判断逆/顺时针
6. 判断线段相交
7. 线段的交点
8. 圆与直线的交点
9. 圆与圆的交点
10. 点的内包
11. 安德鲁算法求凸包
这套算法中, 最基本的是内积和外积, 其他的实现基本是这两个数学概念用法的拓展. 需要进行一些演算.
代码略长, 可自取对应函数(有注释)
PS. 主函数中的数据输入是求凸包.
#include <algorithm>
#include <iostream>
#include <vector>
#include <cmath>
#include <cassert>
#define EPS (1e-10)
#define equals(a, b) (fabs((a) - (b)) < EPS)
using namespace std;
class Point {
public:
double x, y;
Point (double x = 0.0, double y = 0.0): x(x), y(y) {} // 可传递参数的默认构造函数
Point operator + (Point &p) { // 自定义点/向量加法
return Point(x + p.x, y + p.y);
}
Point operator - (Point &p) { // 自定义点/向量减法
return Point(x + p.x, y - p.y);
}
Point operator * (double k) { // 自定义向量乘法(数乘)
return Point(x * k, y * k);
}
Point operator / (double k) {
return Point(x / k, y / k);
}
bool operator < (const Point &p) const {
return x != p.x ? x < p.x : y < p.y;
}
bool operator == (const Point &p) const {
return fabs(x - p.x) < EPS && fabs(y - p.y) < EPS;
}
};
class Segment {
public:
Point p1, p2;
};
typedef Point Vector;
typedef Segment Line;
typedef vector<Point> Polygon;
class Circle {
public:
Point c;
double r;
Circle(Point c = Point(), double r = 0.0): c(c), r(r) {}
};
double norm(Point a) { // 获取向量大小的平方
return a.x * a.x + a.y * a.y;
}
double abs(Point a) { // 获取向量大小, 直接调用norm()
return sqrt(norm(a));
}
double dot(Vector a, Vector b)
{
return a.x * b.x + a.y * b.y;
}
double cross(Vector a, Vector b)
{
return a.x * b.y - a.y * b.x;
}
bool isOrhtogonal(Vector a, Vector b) // 三个重载函数判断是否正交. 利用dot(点积)判断
{
return equals(dot(a, b), 0.0);
}
bool isOrhtogonal(Point a1, Point a2, Point b1, Point b2)
{
return isOrhtogonal(a1 - a2, b1 - b2);
}
bool isOrhtogonal(Segment s1, Segment s2)
{
return isOrhtogonal(s1.p2 - s1.p1, s2.p2 - s2.p1);
}
bool isParallel(Vector a, Vector b) // 三个重载函数判断是否平行.
{
return equals(cross(a, b), 0.0);
}
bool isParallel(Point a1, Point a2, Point b1, Point b2)
{
return isParallel(a1 - a2, b1 - b2);
}
bool isParallel(Segment s1, Segment s2)
{
return isParallel(s1.p2 - s1.p1, s2.p2 - s2.p1);
}
Point Project(Segment s, Point p) // 一个点在一条线段上的投影
{
Vector base = s.p2 - s.p1;
double r = dot(base, p - s.p1) / norm(base);
return base*r + s.p1; // 不能反了, 因为传递的参数是引用; 或者构造一个新的Point变量
}
Point reflect(Segment s, Point p) // 点p关于线段s的对称点, 即映像
{
Point x = Project(s, p);
return (x - p) * 2.0 + p;
}
double getDistance(Point p1, Point p2) // 两点之间的距离
{
return abs(p1 - p2);
}
double getDistanceLP(Line l, Point p) // 点与直线的距离, 直接利用原点和投影点向量
{
return abs(Project(l, p) - p);
}
double getDistanceSP(Segment s, Point p) // 点与线段的距离, 需要做两个特判.
{
if (dot(s.p2 - s.p1, p - s.p1) < 0) return getDistance(s.p1, p);
if (dot(s.p1 - s.p2, p - s.p2) < 0) return getDistance(s.p2, p);
return getDistanceLP(s, p);
}
int ccw(Point p0, Point p1, Point p2) // 判断方向: 逆时针, 顺时针, 射线两端, 线段上.
{
Vector a = p1 - p0;
Vector b = p2 - p0;
if (cross(a, b) > EPS) return 1;
if (cross(a, b) < -EPS) return -1;
if (dot(a, b) < -EPS) return 2;
if (norm(a) < norm(b)) return -2;
return 0;
}
bool intersect(Segment s1, Segment s2) // 利用线段与点之间的方向关系, 判断线段是否相交k
{
return (ccw(s1.p1, s2.p1, s2.p2) * ccw(s1.p2, s2.p1, s2.p2) <= 0) &&
(ccw(s2.p1, s1.p1, s1.p2) * ccw(s2.p2, s1.p1, s1.p2) <= 0);
}
bool intersect(Circle c, Line l) // 相交的重载函数, 判断圆和直线是否相交
{
return abs(Project(l, c.c) - c.c) < c.r;
}
bool intersect(Circle c1, Circle c2) // 判断两圆相交的函数
{
return abs(c1.c - c2.c) < (c1.r + c2.r);
}
double getDistanceSS(Segment s1, Segment s2) // 线段之间的距离, 如果相交则距离为0;
{
if (intersect(s1, s2)) return 0.0;
return min(min(getDistanceSP(s1, s2.p1), getDistanceSP(s1, s2.p2)),
min(getDistanceSP(s2, s1.p1), getDistanceSP(s2, s1.p2)));
}
Point getCrossPoint(Segment s1, Segment s2) // 获取两条线段的交, 前提是这两条线段有交点
{
Vector base = s1.p2 - s1.p1;
double d1 = abs(cross(base, s2.p2 - s1.p1));
double d2 = abs(cross(base, s2.p1 - s1.p1));
double t = d1 / (d1 + d2);
return (s2.p1 - s2.p2) * t + s2.p2;
}
pair<Point, Point> getCrossPoints(Circle c, Line l) // 圆与直线的交点.
{
assert(intersect(c, l)); // 先断言圆和直线是相交的.
Point x = Project(l, c.c);
double d = sqrt(c.r*c.r - norm(x - c.c)) / abs(x - l.p1);
return make_pair((x - l.p1)*d + x, (x - l.p1)*-d + x);
}
double arg(Vector p) { return atan2(p.y, p.x); }
Vector polar(double a, double r) { return Point(cos(r) * a, sin(r) * a); }
pair<Point, Point> getCrossPoints(Circle c1, Circle c2) // 圆与圆之间的两个交点.
{
assert(intersect(c1, c2)); // 断言两圆相交.
double d = abs(c1.c - c2.c);
double a = acos((d * d + c1.r * c1.r - c2.r * c2.r) / (2.0 * d * c1.r));
double t = arg(c2.c - c1.c);
return make_pair(polar(c1.r, a + t), polar(c1.r, t - a));
}
int contains(Polygon g, Point p) // 判断是否内包, ON 1; IN 2; OUT 0;
{
int n = g.size();
bool x = false;
for (int i = 0; i < n; ++i) {
Point a = g[i] - p, b = g[(i + 1)%n] - p;
if (abs(cross(a, b)) < EPS && dot(a, b) < EPS) return 1;
if (a.y > b.y) swap(a, b);
if (a.y < EPS && EPS < b.y && cross(a, b) > EPS) x = !x;
}
return (x ? 2 : 0);
}
bool find(Polygon s, Point p)
{
for (size_t i = 0; i < s.size(); ++i) {
if (s[i] == p) return true;
}
return false;
}
Polygon andrewScan(Polygon g) // 求凸包的安德鲁算法.
{
if (g.size() < 3) return g; // 边界特殊情况
sort(g.begin(), g.end());
int n = g.size();
Polygon s;
s.push_back(g[0]);
Point one = g[0];
s.push_back(g[1]);
Point two = g[1];
for (int i = 2; i < n; ++i) {
if (ccw(one, two, g[i]) == 1) { // 当为逆时针时弹出
s.pop_back();
}
s.push_back(g[i]);
one = s[s.size() - 2];
two = s.back();
}
Polygon v;
v.push_back(g.back());
int k = 0;
for (size_t i = g.size() - 2; i >= 0; --i) {
if (!find(s, g[i])) {
k = i;
v.push_back(g[i]);
break;
}
}
one = v[0];
two = v[1];
for (int i = k - 1; i >= 0; --i) {
if (ccw(one, two, g[i]) == -1) { // 当为逆时针时
v.pop_back();
}
v.push_back(g[i]);
one = v[v.size() - 2];
two = v.back();
}
for (int i = 1; i < (int)v.size() - 1; ++i) {
s.push_back(v[i]);
}
reverse(s.begin() + 1, s.end());
return s;
}
int main()
{
Polygon p;
int n;
cin >> n;
for (int i = 0; i < n; ++i) {
int x, y;
cin >> x >> y;
p.push_back(Point(x, y));
}
Polygon a = andrewScan(p);
cout << a.size() << endl;
for (int i = 0; i < (int)a.size(); ++i) {
cout << a[i].x << ' ' << a[i].y << endl;
}
}
/*
7
2 1
0 0
1 2
2 2
4 2
1 3
3 3
*/