计算几何模版


计算几何:

极角排序:

#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
struct node
{
    double x, y;
} a[55];
//计算叉积p0p1×p0p2
int cross(node p0,node p1,node p2){
    return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
//计算p1p2的距离
double dis(node p1,node p2){
    return sqrt((double)(p2.x-p1.x)*(p2.x-p1.x)+(p2.y-p1.y)*(p2.y-p1.y));
}
//极角排序函数,角度相同则距离小的在前面
bool cmp(node p1,node p2){
    int tmp=cross(a[0],p1,p2);
    if(tmp>0) return true;
    else if(tmp==0&&dis(a[0],p1)<dis(a[0],p2)) return true;
    else return false;
}
int main()
{
    int i = 0;
    while (scanf("%lf%lf", &a[i].x, &a[i].y) != EOF)
    {
        i++;
    }
    //cout << i << endl;
    sort(a + 1, a + i, cmp);
    for (int k = 0; k < i; k++)
    {
        printf("(%d,%d)\n", (int)a[k].x, (int)a[k].y);
    }
    return 0;
}

凸包—Graham算法:

#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstdio>
using namespace std;
const int maxn = 1e5 + 10;
int n, tot;
struct node
{
    double x, y;
} a[maxn], p[maxn];
double dis(node a, node b)
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
double x(node a, node b, node c)
{
    return (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);
}
bool cmp(node A, node B)
{
    double pp = x(A, B, a[1]);
    if (pp > 0)
        return 1;
    if (pp < 0)
        return 0;
    return dis(a[1], A) < dis(a[1], B);
}
void graham()
{
    int k = 0;
    for (int i = 1; i <= n; i++)
    {
        if (a[i].y < a[k].y || (a[i].y == a[k].y && a[i].x < a[k].x))
            k = i;
    }
    swap(a[1], a[k]);
    sort(a + 2, a + n + 1, cmp);
    tot = 2;
    p[0] = a[0];
    p[1] = a[1];
    for (int i = 2; i <= n; i++)
    {
        while (tot > 0 && x(p[tot - 2], p[tot - 1], a[i]) <= 0)
            tot--;
        p[tot++] = a[i];
    }
    p[tot+1] = a[1];
}
int main()
{
    cin >> n;
    for (int i = 1; i <= n; i++)
    {
        cin >> a[i].x >> a[i].y;
    }
    graham();
    double ans = 0;
    for (int i = 1; i <= tot; i++)
    {
        ans += dis(p[i], p[i + 1]);
        cout << ans << endl;
    }
    printf("%.2lf\n",ans);
    return 0;
}

凸包—Andrew算法:

#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#define ll long long
#define INF 2147483647
#define eps 1e-8
using namespace std;

double ans;
ll n,top,st[100100];
struct node{
	double x,y;
}a[100100];

inline bool cmp(node a,node b){
	return a.x==b.x?a.y<b.y:a.x<b.x;
}

inline double getk(ll i,ll j){
	if(fabs(a[j].x-a[i].x)<eps) return INF;//两个点横坐标相同则斜率为无限大
	return (a[j].y-a[i].y)*1.0/(a[j].x-a[i].x);
}

inline double sqr(double v){
	return v*v;
}

inline double calc(ll i,ll j){
	return sqrt(sqr(a[i].x-a[j].x)+sqr(a[i].y-a[j].y));
}

int main(){
	scanf("%lld",&n);
	for(ll i=1; i<=n; i++) scanf("%lf %lf",&a[i].x,&a[i].y);
	sort(a+1,a+1+n,cmp);
	for(ll i=1; i<=n; i++){
		st[++top]=i;
		while(top>=3&&getk(st[top-2],st[top])<getk(st[top-2],st[top-1])-eps){
			st[top-1]=st[top];
			st[top--]=0;
			//栈中第二个元素出栈
		}
	}
	for(ll i=1; i<top; i++) ans+=calc(st[i],st[i+1]);
	while(top) st[top--]=0;
    
	//同理
	for(ll i=n; i>=1; i--){
		st[++top]=i;
		while(top>=3&&getk(st[top-2],st[top])<getk(st[top-2],st[top-1])-eps){//要求:栈中至少还有三个元素并且斜率符合要求
			st[top-1]=st[top];
			st[top--]=0;
		}
	}
	for(ll i=1; i<top; i++) ans+=calc(st[i],st[i+1]);
	printf("%.2lf",ans);
	return 0;
}

离散化扫描线:

#include <iostream>
#include <algorithm>
#include <string.h>
using namespace std;
#define int long long
const int maxn = 1e6 + 10;

int rk[maxn];
int v[maxn];

struct tree
{
    int left, right;
    int cover;
    int len;
} tr[maxn << 3];

struct L
{
    int x;
    int y1, y2;
    int state;
    bool operator<(L other)
    {
        return x < other.x;
    }
} line[maxn];

void pushup(int pos)
{
    if (tr[pos].cover)
        tr[pos].len = v[tr[pos].right + 1] - v[tr[pos].left];
    else
        tr[pos].len = tr[pos * 2].len + tr[pos * 2 + 1].len;
}

void build(int pos, int left, int right)
{
    tr[pos].left = left;
    tr[pos].right = right;
    if (left == right)
        return;
    int mid = left + (right - left) / 2;
    build(pos * 2, left, mid);
    build(pos * 2 + 1, mid + 1, right);
}

void update(int pos, int left, int right, int w)
{
    if (left <= tr[pos].left && tr[pos].right <= right)
    {
        tr[pos].cover += w;
        pushup(pos);
        return;
    }
    int mid = tr[pos].left + (tr[pos].right - tr[pos].left) / 2;
    if (left <= mid)
        update(pos * 2, left, right, w);
    if (mid < right)
        update(pos * 2 + 1, left, right, w);
    pushup(pos);
}

signed main()
{
    int n;
    cin >> n;
    for (int i = 1; i <= n; i++)
    {
        int a, b, c, d;
        cin >> a >> b >> c >> d;
        line[i] = (L){a, b, d, 1};
        line[i + n] = (L){c, b, d, -1};
        rk[i] = b;
        rk[i + n] = d;
    }
    n *= 2;
    sort(rk + 1, rk + 1 + n);
    int size = unique(rk + 1, rk + n + 1) - rk - 1;
    for (int i = 1; i <= n; i++)
    {
        int p1 = lower_bound(rk + 1, rk + 1 + size, line[i].y1) - rk;
        int p2 = lower_bound(rk + 1, rk + 1 + size, line[i].y2) - rk;
        v[p1] = line[i].y1;
        v[p2] = line[i].y2;
        line[i].y1 = p1;
        line[i].y2 = p2;
    }
    sort(line + 1, line + n + 1);
    build(1, 1, n);
    int ans = 0;
    for (int i = 1; i <= n; i++)
    {
        ans += tr[1].len * (line[i].x - line[i - 1].x);
        update(1, line[i].y1, line[i].y2 - 1, line[i].state);
    }
    cout << ans << endl;
    return 0;
}

随机增量法(最小圆形覆盖):

#include <iostream>
#include <algorithm>
#include <random>
#include <cmath>
using namespace std;
const int maxn = 1e6 + 10;
struct no
{
    double x, y;
} node[maxn];
struct circle
{
    double x, y, r;
};
circle get(no p1, no p2, no p3)
{
    circle p;
    double a1 = p2.x - p1.x, b1 = p2.y - p1.y;
    double c1 = (p2.x * p2.x - p1.x * p1.x + p2.y * p2.y - p1.y * p1.y) / 2;
    double a2 = p3.x - p1.x, b2 = p3.y - p1.y;
    double c2 = (p3.x * p3.x - p1.x * p1.x + p3.y * p3.y - p1.y * p1.y) / 2;
    p.x = (c1 * b2 - c2 * b1) / (a1 * b2 - a2 * b1);
    p.y = (a2 * c1 - a1 * c2) / (b1 * a2 - b2 * a1);
    p.r = sqrt((p.x - p1.x) * (p.x - p1.x) + (p.y - p1.y) * (p.y - p1.y));
    return p;
}
double dis(no a, no b)
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
bool cheak(no p, circle a)
{
    if (sqrt((p.x - a.x) * (p.x - a.x) + (p.y - a.y) * (p.y - a.y)) <= a.r)
        return 1;
    else
        return 0;
}
int main()
{
    int n;
    cin >> n;
    for (int i = 1; i <= n; i++)
    {
        cin >> node[i].x >> node[i].y;
    }
    random_shuffle(node + 1, node + n + 1);
    circle p;
    p.r = 0, p.x = node[1].x, p.y = node[1].y;
    for (int i = 1; i <= n; i++)
    {
        if (!cheak(node[i], p))
        {
            p.r = 0;
            p.x = node[i].x, p.y = node[i].y;
            for (int j = 1; j <= i - 1; j++)
            {
                if (!cheak(node[j], p))
                {
                    p.r = (dis(node[i], node[j])) / 2;
                    p.x = (node[i].x + node[j].x) / 2, p.y = (node[i].y + node[j].y) / 2;
                    for (int k = 1; k <= j - 1; k++)
                    {
                        if (!cheak(node[k], p))
                        {
                            p = get(node[i], node[j], node[k]);
                        }
                    }
                }
            }
        }
    }
    printf("%.2lf %.2lf %.2lf\n", p.x, p.y, p.r);
    return 0;
}

旋转卡壳找最远点对:

#include <iostream>
#include <algorithm>

#include <cmath>
using namespace std;
#define int long long
const int maxn=1e6+10;
int n, cnt; // cnt 是凸包上点的数量
struct node
{                     //结构体,存储坐标系中的点
    double x, y;      //坐标
} d[maxn], s[maxn]; // d 为原点集,s 为凸包上的点集
double check(node p1, node p2, node p3, node p4)
{ //求叉积,也用于下面的 cmp 中
    double x1 = p2.x - p1.x;
    double y1 = p2.y - p1.y;
    double x2 = p4.x - p3.x;
    double y2 = p4.y - p3.y;
    return x1 * y2 - x2 * y1;
}
double dis(node l, node r)
{ //距离公式·平方和
    return (r.x - l.x) * (r.x - l.x) + (r.y - l.y) * (r.y - l.y);
}
bool cmp(node l, node r)
{ //按与最左下点的极角排序
    double tmp = check(d[1], l, d[1], r);
    if (tmp > 0)
        return 1;
    if (tmp == 0 && dis(d[1], l) < dis(d[1], r))
        return 1;
    return 0;
}
void graham()
{
    s[1] = d[1];                 //这样处理下来的第一个点(最左下的点)必然在凸壳上
    sort(d + 2, d + n + 1, cmp); //排序
    for (int i = 2; i <= n; i++)
    {
        while (cnt > 1 && check(s[cnt - 1], s[cnt], s[cnt], d[i]) <= 0)
            cnt--; //叉积小于等于 0 时退栈。
        //注意,小于等于 0,因为等于 0 时相当于存在共线,此时我们求直径就要考虑线上最远的点也就是端点,而不是线上的点。
        cnt++;
        s[cnt] = d[i];
    }
    s[++cnt] = d[1]; //最后一个点连接到第一个点上,让凸包形成闭环
    //凸包构建完毕
}

double get_longest()
{
    double ans = 0;
    int p = 2;     //开始根据边枚举凸壳上的点。我们从第二个点开始枚举。
    s[0] = s[cnt]; //转一周会经过第 n 个点,我们接下来会通过取模避免这一点,但这里就要给s[0]赋值为s[cnt]
    for (int i = 0; i <= cnt; i++)
    {
        while (check(s[p], s[i], s[p], s[i + 1]) < check(s[p + 1], s[i], s[p + 1], s[i + 1]))
            p = (p + 1) % cnt;                                     //枚举凸壳上的点,直到下一个点到边的长度比当前点到边的长度小
        ans = max(ans, max(dis(s[i], s[p]), dis(s[i + 1], s[p]))); //更新答案
    }
    return ans;
}
signed main()
{
    cin >> n;
    node Exchange;           //用于实现交换的工具点
    cin >> d[1].x >> d[1].y; //初始化第一个点
    for (int i = 2; i <= n; i++)
    {
        cin >> d[i].x >> d[i].y;
        if (d[i].y < d[1].y)
        {                    //如果当前点的y坐标小于第一个,我们就让它变成第一个
            Exchange = d[i]; //注意上面的判断是 < 而不是 <= ,因为我们还要保证第一个点是最左下的点
            d[i] = d[1];
            d[1] = Exchange;
        }
    }
    graham();
    double ans = get_longest();//旋转卡壳(平面最远点对)
    printf("%.0lf", ans); //输出答案
    return 0;
}

旋转卡壳找最大矩形覆盖

#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
using namespace std;
const double eps = 1e-11;
const double PI = acos(-1);

int dcmp(double x)
{
    if (fabs(x) < eps)
        return 0;
    return x > 0 ? 1 : -1;
}
struct Point
{
    double x, y;
    Point(double a = 0, double b = 0) : x(a), y(b) {}
};
typedef Point Vector;
typedef vector<Point> Polygon;

Vector operator+(const Vector &a, const Vector &b) { return Vector(a.x + b.x, a.y + b.y); }
Vector operator-(const Vector &a, const Vector &b) { return Vector(a.x - b.x, a.y - b.y); }
Vector operator*(const Vector &a, double &b) { return Vector(a.x * b, a.y * b); }
Vector operator/(const Vector &a, double &b) { return Vector(a.x / b, a.y / b); }
bool operator==(const Vector &a, const Vector &b) { return !dcmp(a.x - b.x) && !dcmp(a.y - b.y); }
bool operator<(const Vector &a, const Vector &b) { return a.x < b.x || (a.x == b.x && a.y < b.y); }
double Dot(const Vector &a, const Vector &b) { return a.x * b.x + a.y * b.y; }
double Length(const Vector &a) { return sqrt(Dot(a, a)); }
double Length2(const Vector &a) { return Dot(a, a); }
double Cross(const Vector &a, const Vector &b) { return a.x * b.y - a.y * b.x; }
double Angle(const Vector &a, const Vector &b) { return acos(Dot(a, b) / Length(a) / Length(b)); }
Vector Rotate(Vector A, double rad) { return Point(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad)); }

int ConvexHull(Point *P, int cnt, Point *res)
{
    sort(P, P + cnt);
    cnt = unique(P, P + cnt) - P;
    int m = 0;
    for (int i = 0; i < cnt; i++)
    {
        while (m > 1 && Cross(res[m - 1] - res[m - 2], P[i] - res[m - 2]) <= 0)
            m--;
        res[m++] = P[i];
    }
    int k = m;
    for (int i = cnt - 2; i >= 0; i--)
    {
        while (m > k && Cross(res[m - 1] - res[m - 2], P[i] - res[m - 2]) <= 0)
            m--;
        res[m++] = P[i];
    }
    if (cnt > 1)
        m--;
    return m;
}

void GetMinRCarea(Point *P, int cnt, double &S, double &C)
{
    S = C = 1e18;
    int l = 1, r = 1, u = 1;
    for (int i = 0; i < cnt; i++)
    { //枚举底边P[i]~P[i+1]
        while (Cross(P[i + 1] - P[i], P[u + 1] - P[u]) > eps)
            u = (u + 1) % cnt; //找上边,夹角还能增加就向后找
        while (Dot(P[i + 1] - P[i], P[r + 1] - P[r]) > eps)
            r = (r + 1) % cnt; //找右边,夹角还小于90°就继续找
        if (i == 0)
            l = (r + 1) % cnt;
        while (Dot(P[i + 1] - P[i], P[l + 1] - P[l]) < -eps)
            l = (l + 1) % cnt;                                                                  //找左边,夹角还大于90°就继续找
        double d = Length(P[i + 1] - P[i]);                                                     //底边长度
        double a = Cross(P[i + 1] - P[i], P[u] - P[i]) / d;                                     //叉积几何意义就是求向量组成的平行四边形面积,除以底边就是高
        double b = (Dot(P[r] - P[i], P[i + 1] - P[i]) - Dot(P[l] - P[i], P[i + 1] - P[i])) / d; //向左向右的投影和(点积几何意义:A*B = |A||B|cosθ一个向量在另一个向量方向上的投影长度乘后者的长度),所以最后还要除后者的长度。
        S = min(S, a * b);
        C = min(C, 2.0 * (a + b));
    }
}
const int N = 100005;
Point P[N], R[N];
int n;

int main()
{
    int t;
    scanf("%d",&t);
    while (t--)
    {
        scanf("%d",&n);
        for (int i = 0; i < n; i++)
            scanf("%lf%lf", &R[i].x, &R[i].y);
        int cnt = ConvexHull(R, n, P);
        P[cnt] = P[0];
        double S, C;
        GetMinRCarea(P, cnt, S, C);
        printf("%.2lf %.2lf\n", S, C);
    }
    return 0;
}

PIck定理:

平面最近点对:

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <set>
const int N = 200005;
int n;
double ans = 1e20;

struct point {
  double x, y;

  point(double x = 0, double y = 0) : x(x), y(y) {}
};

struct cmp_x {
  bool operator()(const point &a, const point &b) const {
    return a.x < b.x || (a.x == b.x && a.y < b.y);
  }
};

struct cmp_y {
  bool operator()(const point &a, const point &b) const { return a.y < b.y; }
};

inline void upd_ans(const point &a, const point &b) {
  double dist = sqrt(pow((a.x - b.x), 2) + pow((a.y - b.y), 2));
  if (ans > dist) ans = dist;
}

point a[N];
std::multiset<point, cmp_y> s;

int main() {
  scanf("%d", &n);
  for (int i = 0; i < n; i++) scanf("%lf%lf", &a[i].x, &a[i].y);
  std::sort(a, a + n, cmp_x());
  for (int i = 0, l = 0; i < n; i++) {
    while (l < i && a[i].x - a[l].x >= ans) s.erase(s.find(a[l++]));
    for (auto it = s.lower_bound(point(a[i].x, a[i].y - ans));
         it != s.end() && it->y - a[i].y < ans; it++)
      upd_ans(*it, a[i]);
    s.insert(a[i]);
  }
  printf("%.4lf", ans);
  return 0;
}

平面最小正方形覆盖(三分):

#include<iostream>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
#define INF 0x3f3f3f3f
#define PI acos(-1.0)
using namespace std;
const double eps =1e-12;
const int maxn = 50;
struct Node{
	int x;
	int y;
}a[maxn];
int n;
double C(double degree){
	double x,y,maxx,maxy,minx,miny;
	maxx = maxy = -INF;
	minx = miny = INF;
	for(int i=0;i<n;i++){
		x=a[i].x*cos(degree)-a[i].y*sin(degree);
		y=a[i].x*sin(degree)+a[i].y*cos(degree);
		maxx=max(maxx,x);
		maxy=max(maxy,y);
		minx=min(minx,x);
		miny=min(miny,y);
	}
	return max(maxx-minx,maxy-miny);
}
double bsearch(){
	double left=0,right=PI/2.0;
	double midl,midr;
	while(right-left>eps){
		midl=(left+right)/2;
		midr=(midl+right)/2;
		if(C(midl)<C(midr)){
			right=midr;
		}
		else{
			left=midl;
		}
	}
	return C(left);
}
int main()
{
	int T;
	cin>>T;
	while(T--){
		cin>>n;
		for(int i=0;i<n;i++){
			cin>>a[i].x>>a[i].y;
		}
		double low = bsearch();
		printf("%.2lf\n",low*low);
	}
	return 0;
}

文章作者: fatzard
版权声明: 本博客所有文章除特別声明外,均采用 CC BY 4.0 许可协议。转载请注明来源 fatzard !
评论
  目录
本站总访问量