[好题][半平面交][dfs]Nice Milk POJ1271

Little Tomy likes to cover his bread with some milk. He does this by putting it in the cup so that one of its sides (called bottom side) touches the bottom of the cup, just as the picture below:

Since the milk in the cup is limited, only part of the bread is covered with milk(as shown in the pictures). That is, only the area between the surface of the milk and the bottom side of the bread is covered. Note that the distance between these two lines is always h - the depth of the milk, which is also known to him.
Tomy wants to cover this bread with largest possible area of milk in this way, but he doesn't want to do more than k actions. Help him, will you?
(You may assume that the cup is wide enough, wider than any side of the bread, so it's possible to cover any side completely)

Input

The input will contain no more than 10 test cases. Each test case begins with a line containing three integers n,k and h (3<=n<=20, 0<=k<=8, 0<=h<=10). A piece of bread is guaranteed to be a convex polygon of n vertices. In the following n lines, each line contains two integers xi and yi(0<=xi,yi<=1000), representing the Cartesianism coordinate of the ith vertex. The vertices are anti-clockwise numbered. The test case containing n=0, k=0, h=0 will terminate the input, you should not give an answer to this case.

Output

Output the area of the largest possible area of bread covered with milk with two decimal places. Output one line for each test case.

Sample Input

4 2 1
1 0
3 0
5 2
0 4
0 0 0

Sample Output

7.46

题意： 有一块凸多边形面包，它有n条边，现在你可以把它泡进牛奶k次，牛奶高度为h，询问k次后有多少面积被牛奶浸泡。

分析： 毒瘤题，卡时间卡的很死，稍微不注意细节就超时了。现在总结一下TLE一整天的经验。 首先如果已经知道哪些边要被泡进牛奶然后求面积是很简单的，直接让那些边推进h长度然后求半平面交的面积，最后拿总面积一减就是牛奶浸泡面积。由于不知道哪k条边要被浸泡，很自然想到可以dfs出来所有组合情况，于是随手写了个dfs，很自然地TLE了。

优化部分：首先dfs可以进行优化，最好想的剪枝就是如果到当前位置时后面的边选不到k条那就直接退出。还可以想到如果牛奶浸泡面积已经等于面包面积了，之后也就不需要选边了，答案就是面包面积。dfs部分优化完了，之后是求半平面交的优化，在求各边半平面交时需要先对各边排序+去重，正常情况下每次dfs到终点都需要进行一次求半平面交，也就是会进行巨多次对边的排序+去重，实际上给出面包形状后各边的倾斜角已经固定，只需要最开始排序一次就够了，另外面包是凸多边形，去重根本去不掉任何边，所以不需要去重。还有一个优化是预处理出每条边的法向量，因为在dfs中将各边推进时每次都要用到这个法向量，如果现场求的话太浪费时间，预处理出来每次直接调用就行。最后一个优化是n = k的情况，这时候就不用dfs了，直接所有边都泡就行。

如果wa了可能是没考虑k > n的情况，此时dfs跑不到终点，没法正确更新答案，直接让k = min(n, k)。

具体代码如下：

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
using namespace std;
const double eps = 1e-8;
const double inf = 1e20;
const double pi = acos(-1.0);
const int maxp = 25;
int sgn(double x)
{
	if(fabs(x) < eps)return 0;
	if(x < 0)return -1;
	else return 1;
}
struct Point
{
	double x, y;
	Point(){}
	Point(double _x,double _y){x = _x, y = _y;}
	void input(){scanf("%lf%lf",&x,&y);}
	void output(){printf("%.2f %.2f\n",x,y);}
	bool operator == (Point b)const{return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;}
	bool operator < (Point b)const{return sgn(x-b.x)== 0?sgn(y-b.y)<0:x<b.x;}
	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;}
	//返回长度
	double len(){return sqrt(x*x+y*y);/*库函数*/}
	//返回长度的平方
	double len2(){return x*x + y*y;}
	//返回两点的距离
	double distance(Point p){return hypot(x-p.x,y-p.y);}
	Point operator +(const Point &b)const{return Point(x+b.x,y+b.y);}
	Point operator *(const double &k)const{return Point(x*k,y*k);}
	Point operator /(const double &k)const{return Point(x/k,y/k);}
}; 
struct Line
{
	Point s,e;
	Line(){}
	Line(Point _s,Point _e){s = _s, e = _e;}
	bool operator ==(Line v){return (s == v.s)&&(e == v.e);}
	//`根据一个点和倾斜角angle确定直线,0<=angle<pi`
	Line(Point p,double angle)
	{
		s = p;
		if(sgn(angle-pi/2) == 0){e = (s + Point(0,1));}
		else{e = (s + Point(1,tan(angle)));}
	}
	void input()
	{
		s.input();
		e.input();
	}
	bool parallel(Line v){return sgn((e-s)^(v.e-v.s)) == 0;/*两向量叉积为0*/ }
	Point crosspoint(Line v)
	{
		double a1 = (v.e-v.s)^(s-v.s);
		double a2 = (v.e-v.s)^(e-v.s);
		return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
	}
};
struct halfplane:public Line
{
	double angle;
	halfplane(){}
	//`表示向量s->e逆时针(左侧)的半平面`
	halfplane(Point _s,Point _e)
	{
		s = _s;
		e = _e;
	}
	halfplane(Line v)
	{
		s = v.s;
		e = v.e;
	}
	void calcangle(){angle = atan2(e.y-s.y,e.x-s.x);}
	bool operator <(const halfplane &b)const{return angle < b.angle;}
};
struct polygon
{
	int n;
	Point p[maxp];
	Line l[maxp];
	double getarea()
	{
		double sum = 0;
		for(int i = 0;i < n;i++){
			sum += (p[i]^p[i+1]);
		}
		return fabs(sum)/2;
	}
};
struct halfplanes
{
	int n;//需要输入 
	halfplane hp[maxp];//需要输入，且封闭区域都在向量逆时针方向 
	Point p[maxp];
	int que[maxp];
	int st,ed;//队列的头尾指针，且下标从0开始，指向元素就是头和尾 
	void push(halfplane tmp){hp[n++] = tmp;}
	//去重
	void unique()
	{
		int m = 1;
		for(int i = 1;i < n;i++)
		{
			if(sgn(hp[i].angle-hp[i-1].angle) != 0)
				hp[m++] = hp[i];
			//去除极角相同的情况下，位置在右边(沿向量方向)的边 
			else if(sgn( (hp[m-1].e-hp[m-1].s)^(hp[i].s-hp[m-1].s) ) > 0)
				hp[m-1] = hp[i];
		}
		n = m;
	}
	bool halfplaneinsert()//判断半平面交是否存在 
	{
		que[st=0] = 0;
		que[ed=1] = 1;
		p[1] = hp[0].crosspoint(hp[1]);
		for(int i = 2;i < n;i++){
			while(st<ed && sgn((hp[i].e-hp[i].s)^(p[ed]-hp[i].s))<0)ed--;
			while(st<ed && sgn((hp[i].e-hp[i].s)^(p[st+1]-hp[i].s))<0)st++;
			que[++ed] = i;
			if(hp[i].parallel(hp[que[ed-1]]))return false;
			p[ed]=hp[i].crosspoint(hp[que[ed-1]]);
		}
		while(st<ed && sgn((hp[que[st]].e-hp[que[st]].s)^(p[ed]-hp[que[st]].s))<0)ed--;
		while(st<ed && sgn((hp[que[ed]].e-hp[que[ed]].s)^(p[st+1]-hp[que[ed]].s))<0)st++;
		if(st+1>=ed)return false;//最后剩下小于三条直线，表明半平面交不存在 
		return true;
	}
	void getconvex(polygon &con)
	{
		p[st] = hp[que[st]].crosspoint(hp[que[ed]]);
		con.n = ed-st+1;
		for(int j = st,i = 0;j <= ed;i++,j++)
			con.p[i] = p[j];
		con.p[con.n] = con.p[0];
	}
};

int n, k, h;
Point p[maxp];
halfplanes a, cpy;
double ans;
int st[maxp], top;
Point vec[maxp];//预处理法向量
 
void dfs(int now, int num)//到当前位置处已经找了num个位置 
{
	if(k-num > n-now || sgn(ans) == 0)
		return;//剪枝 
	if(num == k)
	{
		for(int i = 1; i <= top; i++)
		{
			int t = st[i];
			cpy.hp[t].s = cpy.hp[t].s+vec[t]*h;
			cpy.hp[t].e = cpy.hp[t].e+vec[t]*h;
		}
		if(cpy.halfplaneinsert())
		{
			polygon b;
			cpy.getconvex(b);
			ans = min(ans, b.getarea());
		}
		else
			ans = 0.0;
		cpy = a;
		return; 
	}
	dfs(now+1, num);
	st[++top] = now;
	dfs(now+1, num+1);
	top--;
}

signed main()
{
	while(~scanf("%d%d%d", &n, &k, &h))
	{
		if(n == 0 && k == 0 && h == 0)
			break;
		ans = 1e20;
		for(int i = 0; i < n; i++)
			p[i].input();
		p[n] = p[0];
		a.n = 0;
		for(int i = 0; i < n; i++)
			a.push(halfplane(p[i], p[i+1]));
		for(int i = 0;i < n;i++)
			a.hp[i].calcangle();
		sort(a.hp,a.hp+n);//先对倾斜角排序 
		cpy = a;
		k = min(n, k);//很重要的处理
		for(int i = 0; i < n; i++)//预处理每边的法向量 
		{
			vec[i] = a.hp[i].e-a.hp[i].s;
			vec[i] = Point(-vec[i].y, vec[i].x);
			vec[i] = vec[i]/vec[i].len();
		}
		if(k < n)//这时候才需要dfs选边 
		{
			top = 0;
			dfs(0, 0);
		} 
		else
		{
			for(int i = 0; i < n; i++)
			{
				cpy.hp[i].s = cpy.hp[i].s+vec[i]*h;
				cpy.hp[i].e = cpy.hp[i].e+vec[i]*h;
			}
			cpy.halfplaneinsert();
			polygon con;
			cpy.getconvex(con);
			ans = con.getarea(); 
		}
		a.halfplaneinsert();
		polygon con;
		a.getconvex(con);
		printf("%.2f\n", con.getarea()-ans);
	}
    return 0;
}