POJ 3525 二分答案，推进多边形和半平面交

题目在问这样一个问题：给定一个凸多边形，找到其中的一个点，使得其到每条边的距离最小值最大，输出这个距离。

其实就是在问你，这个多边形中最大的一个内切圆有多大。

怎么做呢？如果我们事先知道一个半径R，我们是不是能验证这个R是否可行呢？

答案是肯定的，这样想：如果我们把这个多边形每条边都向内推进R，之后如果这个多边形还存在的话，就说明这个半径为R的圆肯定塞得下，因为还可以往里缩嘛。直到他们缩成一个点的时候，这个时候缩的距离就是我们要找的R了，同时这个点也肯定是我们要找的圆的圆心。

那么，我们就可以二分一下R，之后用半平面交判断一下多边形是否还在就行了。

半平面交的实现和求多边形的核一模一样，我们可以考虑这个问题是一堆半平面的交，实际上实现手段是切割已有的答案区域，得到一个更小的答案区域，这就是半平面交的实现手段。为此，我们一开始选择一个肯定要大一些的答案区域，再不断用直线去切割就可以了。

最想说一下什么呢？是“将每条边推进R”应该怎么实现。很多人绕了半天，实际上我找到一个很简单的公式：把直线表示成Ax + By + C = 0后（已考虑顺逆时针所带来的正负号问题），把直线向内推进d时，只需要写：

c += d * sqrt(a*a + b*b);

就可以了。

为什么呢？我们可以简单推导一下：

①当直线不垂直于X轴的时候，不妨假设直线向左上方移动（直线方向为右上）：

有如下关系成立：

tanθ = k = - A / B

sinα = -B / sqrt( a*a + b*b)

cosα = A / sqrt(a*a + b*b)

dx = d*cosα

dy = d*sinα

新直线方程为A（x+dx） + B（y+dy） + C = 0

拆开整理，会发现分子上面有一个A² + B² ，和下面的sqrt（A² + B²） 约掉一份，得出△C = d* sqrt（A² + B²）。

②当直线垂直于X轴的时候，解析式为Ax + C = 0， 不妨假设直线向左移动（其方向是向上的），则新直线为A（x + d）+ C =0, 于是有△C = Ad = d * sqrt（a*a+b*b）。

所以，以后遇到推进类问题的时候，只需要先求出ABC，再修改一个C，就可以得到推进之后的直线了。

#include<iostream>
#include<string.h>
#include<math.h>
#include<stdio.h>
#include<algorithm>
#define D(x) cout<<#x<<" "<<x<<endl
using namespace std;

const double EPS = 1e-8;
const double INF = 1e8;
int sgn(double x)
{
	if(fabs(x) < EPS)
		return 0;
	if(x>0)
		return 1;
	return -1;
}

struct Point
{
	double x,y;
	Point(){}
	Point(double _x,double _y)
	{
		x = _x;y = _y;
	}
	Point operator -(const Point &b)const
	{
		return Point(x - b.x,y - b.y);
	}
	//叉积
	double operator ^(const Point &b)const
	{
		return x*b.y - y*b.x;
	}
	//点积
	double operator *(const Point &b)const
	{
		return x*b.x + y*b.y;
	}
	//绕原点旋转角度B（弧度值），后x,y的变化
	void transXY(double B)
	{
		double tx = x,ty = y;
		x = tx*cos(B) - ty*sin(B);
		y = tx*sin(B) + ty*cos(B);
	}
};

Point ori[150];
Point kernel[150];
int p_kernel = 0;

Point tp[150];

int n;

void getLine(Point p1,Point p2,double &a,double &b,double &c)
{
	a = p2.y - p1.y;
	b = p1.x - p2.x;
	c = p2.x*p1.y - p1.x*p2.y;
	return;
}

Point Intersection(Point p1,Point p2,double a,double b,double c)
{
	double u = fabs(a*p1.x + b*p1.y + c);
	double v = fabs(a*p2.x + b*p2.y + c);
	Point t;
	t.x = (p1.x*v + p2.x*u)/(u+v);
	t.y = (p1.y*v + p2.y*u)/(u+v);
	return t;
}

void cut(double a,double b,double c,Point p[],int &cnt)
{
	int tmp = 0;
	for(int i = 1;i <= cnt;i++)
	{
		//当前点在左侧，逆时针的点
		if(a*p[i].x + b*p[i].y + c < EPS)
			tp[++tmp] = p[i];
	else
	{
		if(a*p[i-1].x + b*p[i-1].y + c < -EPS)
			tp[++tmp] = Intersection(p[i-1],p[i],a,b,c);
		if(a*p[i+1].x + b*p[i+1].y + c < -EPS)
			tp[++tmp] = Intersection(p[i],p[i+1],a,b,c);
		}
	}
	
	for(int i = 1;i <= tmp;i++)
		p[i] = tp[i];
	p[0] = p[tmp];
	p[tmp+1] = p[1];
	cnt = tmp;
	return;
}

void init()
{
	memset(ori, 0, sizeof(ori));
	memset(kernel, 0, sizeof(kernel));
	p_kernel = 0;
	return;
}

int main()
{
	while(true)
	{
		scanf("%d", &n);
		if(n==0)
			break;
		
		init();
		int i;
		for(i=1; i<=n;i++)
		{
			scanf("%lf %lf", &ori[i].x, &ori[i].y);
		}
		
		ori[0] = ori[n];
		ori[n+1] = ori[1];
		
		double low = 0, high = INF;
		while(fabs(high - low) > 1e-5)
		{
			double mid = (high+low) /2;
				
			for(i=0; i<=n+1;i++)
				kernel[i] = ori[i];
			p_kernel = n;
			
			for(i=1; i<=n;i++)
			{
				double a,b,c;
				getLine(ori[i], ori[i+1], a, b, c);
				
				c += mid * sqrt(a*a + b*b);
				cut(a,b,c, kernel, p_kernel);
			}
			
		//	D(p_kernel);
			
			if(p_kernel > 0)
				low = mid;
			else
				high = mid;
		}
		
		printf("%lf\n", high);
	}
	//system("pause");
	return 0;
}