HDU2823 The widest road

最新推荐文章于 2018-07-27 10:27:22 发布

qhpeklh5959

最新推荐文章于 2018-07-27 10:27:22 发布

阅读量803

点赞数

CC 4.0 BY-SA版权

分类专栏：计算几何

本文链接：https://blog.youkuaiyun.com/severus_qin/article/details/12350061

计算几何专栏收录该内容

9 篇文章

订阅专栏

本文介绍了一种解决旋转卡壳求凸包最近距离问题的算法实现过程，包括几何运算、凸包计算及优化技巧。通过旋转卡壳板子实例，详细阐述了算法的具体步骤和关键细节。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

点击打开链接

旋转卡壳求凸包最近距离，可能是没判断凸包交点什么的，一直RE……先留着当旋转卡壳板子……

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <vector>
#include <iostream>
using namespace std;
const double eps = 1e-8;
int dblcmp(double x){
	if (fabs(x) < eps) return 0;
	return x > 0 ? 1 : -1;
}
const double pi = acos(-1.0);
inline double sqr(double x){
	return x * x;
}
struct Point2D{
	double x, y;
	Point2D(){}
	Point2D(double _x, double _y) : x(_x), y(_y){}
	void input(){
		scanf("%lf%lf", &x, &y);
	}
	void output(){
		printf("%.2lf %.2lf\n", x, y);
	}
	friend Point2D operator + (const Point2D &a, const Point2D &b){
		return Point2D(a.x + b.x, a.y + b.y);
	}
	friend Point2D operator - (const Point2D &a, const Point2D &b){
		return Point2D(a.x - b.x, a.y - b.y);
	}
	friend Point2D operator * (const double &a, const Point2D &b){
		return Point2D(a * b.x, a * b.y);
	}
	friend Point2D operator / (const Point2D &a, const double &b){
		return Point2D(a.x / b, a.y / b);
	}
	friend bool operator == (const Point2D &a, const Point2D &b){
		return dblcmp(a.x - b.x) == 0 && dblcmp(a.y - b.y) == 0;
	}
	double len(){
		return sqrt(sqr(x) + sqr(y));
	}
};
double det(const Point2D &a, const Point2D &b){
	return a.x * b.y - a.y * b.x;
}
double dot(const Point2D &a, const Point2D &b){
	return a.x * b.x + a.y * b.y;
}
double dist(const Point2D &a, const Point2D &b){
	return (a - b).len();
}
Point2D rotate_point(const Point2D &p, double A){
	double tx = p.x, ty = p.y;
	return Point2D(tx * cos(A) - ty * sin(A), tx * sin(A) + ty * cos(A));
}
const int maxp = 1000 + 10;
int n, m;
bool cmpx(Point2D a, Point2D b){
	return dblcmp(a.x - b.x) < 0 || (dblcmp(a.x - b.x) == 0 && dblcmp(a.y - b.y) < 0);
}
struct Polygon2D{
	int n;
	Point2D p[maxp];
	Polygon2D(){}
	void copy(Polygon2D &A){
		A.n = n;
		for (int i = 0; i < n; i++)
			A.p[i] = p[i];
	}
	void input(){
		for (int i = 0; i < n; i++)
			p[i].input();
	}
	void output(){
		for (int i = 0; i < n; i++)
			p[i].output();
		puts("");
	}
	struct cmp{
		Point2D p;
		cmp(const Point2D &p0){p = p0;}
		bool operator ()(const Point2D &aa, const Point2D &bb){
			Point2D a = aa, b = bb;
			int d = dblcmp(det(a - p, b - p));
			if (d == 0)
				return dblcmp(dist(a, p) - dist(b, p)) < 0;
			return d > 0;
		}
	};
	void norm(){
		Point2D mi = p[0];
		for (int i = 0; i < n; i++)
			if (cmpx(p[i], mi)) mi = p[i];
		sort(p, p + n, cmp(mi));
	}
	void getconvex(Polygon2D &convex){
		sort(p, p + n, cmpx);
		convex.n = n;
		for (int i = 0; i < min(n, 2); i++)
			convex.p[i] = p[i];
		if (n <= 2) return;
		int &top = convex.n;
		top = 1;
		for (int i = 2; i < n; i++){
			while(top && dblcmp(det(convex.p[top] - p[i], convex.p[top - 1] - p[i])) <= 0)
				top -= 1;
			convex.p[++top] = p[i];
		}
		int tmp = top;
		convex.p[++top] = p[n - 2];
		for (int i = n - 3; i >= 0; i--){
			while(top != tmp && dblcmp(det(convex.p[top] - p[i], convex.p[top - 1] - p[i])) <= 0)
				top -= 1;
			convex.p[++top] = p[i];
		}
	}
};
struct Segment2D{
	Point2D p1, p2;
	int sz;
	Segment2D(){sz = 0;}
	Segment2D(Point2D _p1) : p1(_p1){sz = 1;}
	Segment2D(Point2D _p1, Point2D _p2) : p1(_p1), p2(_p2){sz = 2;}
};
struct event{
	double argu;
	Segment2D seg[2];
	event(){}
	event(double _argu, Segment2D _s1, Segment2D _s2) : argu(_argu){
		seg[0] = _s1; seg[1] = _s2;
	}
	bool operator < (const event &A) const {
		return argu < A.argu;
	}
};
double fix(double v){
	return v < 0 ? v + pi : v;
}
double ang(Point2D p){
	if (p.x == 0) return pi / 2;
	return atan2(p.y, p.x);
}
void rot(Polygon2D cvx, vector<event>& evt, vector<double>& evt_arg){
	int n = cvx.n;
	evt.clear();
	int si = 0, sj = 0;
	for (int i = 0; i < n; i++){
		if (!cmpx(cvx.p[i], cvx.p[si])) si = i;
		if (cmpx(cvx.p[i], cvx.p[sj])) sj = i;
	}
	int i = si, j = sj;
	while(i != sj || j != si){
		int vi = i, vj = j;
		double argu;
		Segment2D seg1, seg2;
		double crs = det(cvx.p[(i + 1) % n] - cvx.p[i], cvx.p[(j + 1) % n] - cvx.p[j]);
		if (dblcmp(crs) < 0){
			vi = (i + 1) % n;
			argu = fix(ang(cvx.p[vi] - cvx.p[i]));
			seg1 = Segment2D(cvx.p[vi], cvx.p[i]);
			seg2 = Segment2D(cvx.p[j]);
		}else if (dblcmp(crs) == 0){
			vi = (i + 1) % n;
			vj = (j + 1) % n;
			argu = fix(ang(cvx.p[vi] - cvx.p[i]));
			seg1 = Segment2D(cvx.p[vi], cvx.p[i]);
			seg2 = Segment2D(cvx.p[vj], cvx.p[j]);
		}else{
			vj = (j + 1) % n;
			argu = fix(ang(cvx.p[vj] - cvx.p[j]));
			seg1 = Segment2D(cvx.p[vj], cvx.p[j]);
			seg2 = Segment2D(cvx.p[i]);
		}
		i = vi; j = vj;
		evt.push_back(event(argu, seg1, seg2));
		evt_arg.push_back(argu);
	}
	sort(evt.begin(), evt.end());
}
double point_seg_dist(Point2D a, Point2D b, Point2D c){
	double ans = min(dist(a, b), dist(a, c));
	if (dblcmp(dot(a - c, b - c)) >= 0 && dblcmp(dot(c - b, a - b)) >= 0){
		double val = fabs(det(b - a, c - a) / (b - c).len());
		ans = min(ans, val);
	}
	return ans;
}
double min_dist(Segment2D s1, Segment2D s2){
	if (s1.sz > s2.sz) swap(s1, s2);
	if (s1.sz == 1){
		if (s2.sz == 1) return dist(s1.p1, s2.p1);
		else return point_seg_dist(s1.p1, s2.p1, s2.p2);
	}else{
		double ans = 1e20;
		for (int o = 0; o < 2; o++){
			double dis1 = point_seg_dist(s1.p1, s2.p1, s2.p2);
			double dis2 = point_seg_dist(s1.p2, s2.p1, s2.p2);
			ans = min(ans, min(dis1, dis2));
			swap(s1, s2);
		}
		return ans;
	}
}
event find(vector<event> &evt, double ag){
	int l = 0, r = evt.size();
	int len = r;
	while(l < r){
		int mid = l + r >> 1;
		if (dblcmp(evt[mid].argu - ag) <= 0) l = mid + 1;
		else r = mid;
	}
	l = (l - 1 + len) % len;
	return evt[l];
}
void solve(){
	Polygon2D p, q, pc, qc;
	p.n = n; p.input();
	q.n = m; q.input();
	p.getconvex(pc); q.getconvex(qc);
	pc.norm(); qc.norm();
	vector<event> evta, evtb;
	vector<double> evt;
	evt.clear();
	rot(pc, evta, evt);
	rot(qc, evtb, evt);
	sort(evt.begin(), evt.end());
	double ans = 1e20;
	for (int i = 0; i < evt.size(); i++){
		double argu = evt[i];
		event eve1 = find(evta, argu);
		event eve2 = find(evtb, argu);
		for (int pp = 0; pp < 2; pp++)
			for (int qq = 0; qq < 2; qq++){
				double dis = min_dist(eve1.seg[pp], eve2.seg[qq]);
				if (dis < ans) ans = dis;
			}
	}
	printf("%.4lf\n", ans);
}
int main(){
	#ifndef ONLINE_JUDGE
	freopen("in", "rt", stdin);
	#endif
	while(scanf("%d%d", &n, &m) == 2) solve();
	return 0;
}