[计算几何][SHOI2008]安全的航线flight

Description

在设计航线的时候，安全是一个很重要的问题。首先，最重要的是应采取一切措施确保飞行不会发生任何事故，但同时也需要做好最坏的打算，一旦事故发生，就要确保乘客有尽量高的生还几率。当飞机迫降到海上的时候，最近的陆地就是一个关键的因素。航线中最危险的地方就是距离最近的陆地最远的地方，我们称这种点为这条航线“孤地点”。孤地点到最近陆地的距离被称为“孤地距离”。作为航空公司的高级顾问，你接受的第一个任务就是尽量找出一条航线的孤地点，并计算这条航线的孤地距离。为了简化问题，我们认为地图是一个二维平面，陆地可以用多边形近似，飞行线路为一条折线。航线的起点和终点都在陆地上，但中间的转折点是可能在海上（如下图所示，方格标示出了孤地点）。 

Input

输入的第一行包括两个整数C和N（1≤C≤20，2≤N≤20），分别代表陆地的数目的航线的转折点的数目。接下来有N行，每行有两个整数x,y。(x,y)表示一个航线转折点的坐标，第一个转折点为航线的起点，最后一个转折点为航线的终点。接下来的输入将用来描述C块大陆。每块输入由一个正整数M开始（M≤30），M表示多边形的顶点个数，接下来的M行，每行会包含两个整数x,y，(x,y)表示多边形的一个顶点坐标，我们保证这些顶点以顺时针或逆时针给出了该多边形的闭包，不会出现某些边相交的情况。此外我们也保证输入数据中任何两块大陆不会相交。输入的所有坐标将保证在-10000到10000的范围之间。

Output

输出一个浮点数，表示航线的孤地距离，数据保留2位小数。

二分的算法是很好出的：每次把陆地向外推出一个距离，然后判断所有的航线是否被覆盖，这样做就等价于题目中所描述的孤地点了。

但是，我苦逼地写了 3 遍才 A 。。。。计算几何的代码量真是动辄一百多行啊。。。

先随便乱写了一个，结果写的稀乱而且各种 WA 。。。

然后，我发现多边形的顶点向外推将形成一个圆弧。一开始，我试图直接将多边形用圆弧和线段描述出来，但是发现这样有一个严重的问题：多边形是任意多边形而不是凸多边形。这样一来如果有角大于平角则必然会将两条边重叠地推出，然后圆弧也就乱了，于是直接求交便是不可能的了。

之后我干脆把多边形剖分成原多边形和多个矩形和多个扇形。然后由于是完全覆盖的，所以扇形当作圆处理即可。结果不但好写而且一下就 A 了。。。

错点：

其一，当一条直线穿过多边形的一个顶点时，如果用两条边去交或者不交，都会导致出现奇数个点。。。所以需要用点交或者判重。。。（我是用点交的。。。）

其二，当判断一个点是否在一个线段上，且如不在则按情况分别返回两端点。这里各种平行于坐标轴的线段简直虐翻我了。。。

Code :

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
using namespace std;
#define eps 0.000000001
#define oo 999999999999.

struct point {double x, y;};
struct seg {point x, y;};
struct line {point f; double c;};

int sig(double a) {return a < - eps ? - 1 : a > eps;}
double len(point a) {return sqrt(a.x * a.x + a.y * a.y);}

double operator *(point a, point b) {return a.x * b.y - a.y * b.x;}
double operator %(point a, point b) {return a.x * b.x + a.y * b.y;}

point operator +(point a, point b) {return (point) {a.x + b.x, a.y + b.y};}
point operator -(point a, point b) {return (point) {a.x - b.x, a.y - b.y};}

point operator *(line a, line b)
{
  double c = a.f * b.f;
	return (point) {(b.f.y * a.c - a.f.y * b.c) / c, 
    (a.f.x * b.c - b.f.x * a.c) / c};
}

point make_point(point a, double dist)
{
  double x = a.x * a.x;
  double y = a.y * a.y, z = x + y;
  return (point) {sqrt(x * dist * dist / z) * sig(a.x),
    sqrt(y * dist * dist / z) * sig(a.y)};
}

line make_line(seg a)
{
  point f = (point) {a.y.y - a.x.y, a.x.x - a.y.x};
  return (line) {f, f.x * a.x.x + f.y * a.x.y};
}

point select(seg a, point b)
{
	if (sig(a.x.x - a.y.x))
	{
		if ((a.x.x - b.x) * (a.y.x - b.x) < eps) return b;
		if (a.x.x > a.y.x) swap(a.x, a.y);
		return b.x < a.x.x ? a.x : a.y;
	}
	else
	{
 		if ((a.x.y - b.y) * (a.y.y - b.y) < eps) return b;
		if (a.x.y > a.y.y) swap(a.x, a.y);
		return b.y < a.x.y ? a.x : a.y;
	}
}

point operator *(seg a, line b)
{
  point p = make_line(a) * b;
  if (sig(a.x.x - a.y.x) && (a.x.x - p.x) * (a.y.x - p.x) > - eps) return (point) {oo};
  if (sig(a.x.y - a.y.y) && (a.x.y - p.y) * (a.y.y - p.y) > - eps) return (point) {oo};
  return p;
}

point cross_circle(line a, seg b)
{
	line l = make_line(b);
	line r = (line) {(point) {l.f.y, - l.f.x}, a.f.x * l.f.y - a.f.y * l.f.x};
	point p = l * r, q;
	double dist = len(p - a.f), A, B;
	if (dist - a.c > - eps) return (point) {oo};
	dist = sqrt(a.c * a.c - dist * dist);
	q = make_point(b.x - b.y, dist);
	A = len(select(b, p + q) - b.x), B = len(select(b, p - q) - b.x);
	return A < B ? (point) {A, B} : (point) {B, A};
}

void cross_poly(int k, point poly[], seg b, int & tot, point u[])
{
	int i, j; double a[31];
	seg s; point p; line l = make_line(b); j = 0;
	for (i = 1; i <= k; ++ i)
	{
		s = (seg) {poly[i - 1], poly[i]}, p = s * l;
		if (p.x != oo) a[++ j] = len(select(b, p) - b.x);
		if (! sig((b.x - poly[i]) * (b.y - poly[i])))
			a[++ j] = len(select(b, poly[i]) - b.x);
	}
	sort(a + 1, a + j + 1);
	for (i = 1; i <= j; i += 2)
		u[++ tot] = (point) {a[i], a[i + 1]};
}

int C, n;
int num[35], total[35];
point fli[35];
point lan[35][35], cover[35][1000];

bool cmp(point A, point B) {return A.x < B.x;}

bool okay(double dist)
{
	int i, j, k, tot; double le;
	point p; point cov[1000], poly[6]; seg s;
	for (i = 2; i <= n; ++ i)
	{
		tot = total[i], memcpy(cov, cover[i], sizeof cover[i]);
		for (j = 1; j <= C; ++ j)
		{
			s = (seg) {fli[i - 1], fli[i]};
			for (k = 1; k <= num[j]; ++ k)
			{
				p = cross_circle((line) {lan[j][k], dist}, s);
				if (p.x != oo) cov[++ tot] = p;
				
				p = lan[j][k] - lan[j][k - 1];
				p = make_point((point) {p.y, - p.x}, dist);
				poly[1] = lan[j][k - 1], poly[2] = lan[j][k - 1] + p;
				poly[3] = lan[j][k] + p, poly[0] = poly[4] = lan[j][k];
				cross_poly(4, poly, s, tot, cov);
			}
		}
		sort(cov + 1, cov + tot + 1, cmp);
		p = cov[1], le = 0;
		for (j = 2; j <= tot; ++ j)
			if (cov[j].x > p.y) le += p.y - p.x, p = cov[j];
			else if (cov[j].y > p.y) p.y = cov[j].y;
		le += p.y - p.x;
		if (len(s.y - s.x) - le > eps) return 0;
	}
	return 1;
}

int main()
{
	freopen("1020.in", "r", stdin);
	freopen("1020.out", "w", stdout);
	
	int i, j; double l, r, s;
	scanf("%d%d", & C, & n);
	for (i = 1; i <= n; ++ i)
		scanf("%lf%lf", & fli[i].x, & fli[i].y);
	for (i = 1; i <= C; ++ i)
	{
		scanf("%d", num + i);
		for (j = 1; j <= num[i]; ++ j)
			scanf("%lf%lf", & lan[i][j].x, & lan[i][j].y);
		
		s = lan[i][num[i]] * lan[i][1];
		for (j = 2; j <= num[i]; ++ j)
			s += lan[i][j - 1] * lan[i][j];
		if (s < - eps)
			for (j = 1; j << 1 <= num[i]; ++ j)
				swap(lan[i][j], lan[i][num[i] - j + 1]);
		
		lan[i][0] = lan[i][num[i]];
		for (j = 2; j <= n; ++ j)
			cross_poly(num[i], lan[i], (seg) {fli[j - 1], fli[j]}, total[j], cover[j]);
	}
	for (l = 0, r = 10000; r - l > 0.001; )
		okay(s = (l + r) / 2) ? r = s : l = s;
	printf("%.2lf", s);
	
	return 0;
}