POJ 3525 Most Distant Point from the Sea 二分+内推+半平面交

题目描述：http://poj.org/problem?id=3525 

解题思路：

题目要求的是多边形内的点到边的距离的最小值的最大值。一般看到最小值的最大值和最大值的最小值这种问题第一反应就是二分.. 将最优性问题转化为可行性问题。所以我们可以二分点到边的最小距离D，如果存在点p使得Dp=D，那么可能有另外一个点p'使Dp'>D。

那么我们的问题转化成了，判断多边形内有没有一点使得这个点到边的最小距离为D。

如果这时候在草稿纸上写写画画，很容易就可以看出如果我们把每条边向多边形内推进D，求每条边决定的半平面的交，区域中的点到所有边的距离都>=D，区域边界上的点到某条边距离为D，到其他边距离>=D。所以如果交区域存在则证明有p满足条件，否则没有。

那么问题转化成了，二分最大距离，内推边，求半平面交。

至于内推边，假设有向边端点为(x1, y1)(x2, y2)，令v=(x2-x1, y2-y1)，移动向量m为(x, y)。因为m在v顺时针90°方向，所以m×v=|v|*d；因为v、m夹角为90°，所以v·m=0。联立两个方程解出(x, y)，那么平移后的端点就是(x1+x, y1+y)、(x2+x, y2+y)。

复杂度：

假设二分的范围为0到L，求半平面交为O(n^2)，内推边为O(n)，总复杂度为O(n^2*logL)。

/* 
 *  <span style="white-space:pre">	</span>2014.11.18
 *	Problem: 3525		
 *	Memory: 152K		Time: 0MS
 *	Language: C++		Result: Accepted
 *
 */
#include "iostream"
#include "cmath"
#define EPS 1e-8
#define MAXN 107
#define INF 1e6

struct Point {
	double x, y;
	Point(){}
	Point(double _x, double _y):x(_x),y(_y){}
}p[MAXN];
struct Line {
	Point a, b;
};

Point operator + (const Point &a, const Point &b) {
	return Point(a.x+b.x, a.y+b.y);
}
Point operator - (const Point &a, const Point &b) {
	return Point(a.x-b.x, a.y-b.y);
}

int n, cn;

void init() {
	double x, y;
	for (int i=1; i<=n; i++) scanf("%lf %lf", &p[i].x, &p[i].y);
	p[0] = p[n]; p[n+1] = p[1];
}

void trans(Line t[], double d) {
	Point v, vt;
	double l;
	for (int i=1; i<=n; i++) {
		v = p[i] - p[i-1];
		l = sqrt(pow(v.x, 2) + pow(v.y, 2));
		vt = Point(-v.y*d/l, v.x*d/l);
		t[i].b = p[i] + vt;
		t[i].a = p[i-1] + vt;
	}
	t[0] = t[n]; t[n+1] = t[1];
}

double cross(Point a, Point b, Point c) {  // 叉积
	Point v1 = b-a, v2 = c-a;
	return v1.x*v2.y - v2.x*v1.y;
}

Point intersect(Point p1, Point p2, Point p3, Point p4) { // 直线交点
	double a = p2.y-p1.y, b = -(p2.x-p1.x), c = p1.x*p2.y-p2.x*p1.y;
	double d = p4.y-p3.y, e = -(p4.x-p3.x), f = p3.x*p4.y-p4.x*p3.y;
	return Point((c*e-b*f)/(a*e-b*d), (c*d-a*f)/(b*d-a*e));
}

bool cut(Point c[], Point a, Point b) {  
	Point t[MAXN];
	int tn = 0;
	for (int i=1; i<=cn; i++) {
		if (cross(a, b, c[i])>=-EPS) t[++tn] = c[i];
		else {
			if (cross(a, b, c[i-1])>=EPS) t[++tn] = intersect(a, b, c[i], c[i-1]);
			if (cross(a, b, c[i+1])>=EPS) t[++tn] = intersect(a, b, c[i], c[i+1]);
		}
	}
	if (!tn) return false;
	for (int i=1; i<=tn; i++) c[i] = t[i];
	c[0] = c[cn=tn];
	c[cn+1] = c[1];
	return true;
}

bool planeIntersect(Line p[]) {
	Point c[MAXN];
	c[1] = Point(-INF, -INF); 
	c[2] = Point(-INF, INF);
	c[3] = Point(INF, INF);
	c[4] = Point(INF, -INF);
	cn = 4; c[0] = c[cn]; c[cn+1] = c[1];
	for (int i=1; i<=n; i++) {	// 枚举边， 切割现有凸包点集c, O(n^2)半平面交模板
		if (!cut(c, p[i].a, p[i].b)) return false;
	}
	return true;
}

int main() {
	scanf("%d", &n);
	while (n) {
		init();

		double lans = 0, rans = INF, mid;
		struct Line pt[MAXN];
		// 二分内推的距离
		while (fabs(lans-rans)>=1e-7) {
			mid = (lans+rans)/2;
			trans(pt, mid); 			   // 计算内推后的直线pt[]
			bool fit = planeIntersect(pt); // 检测核是否存在
			if (fit) lans = mid; else rans = mid;
		}
		mid = (lans+rans)/2;
		printf("%lf\n", mid);
		scanf("%d", &n);
	}
	return 0;
}