[AIZU]AIZU - 1146 POJ - 3011 Secrets in Shadows 计算几何

传送门：Secrets in Shadows

题意：
给定\(n(n\le 100)\)个半径为\(1\)的圆柱体，假设太阳光是从无限远处来的平行于地面直线，圆柱会在地面上投射出无限长的矩形阴影。
宽度定义为这些矩形的宽度并，问宽度最大的太阳角度。

题解：
一开始以为是枚举两个圆心，然后垂直连线啥的，发现错了，但是网上找不到好的题解。。于是看std理解了思路，遂写此题解。

显然，我们有了太阳的角度，求宽度是很简单的，只需要逆时针旋转\(\theta\)度，然后计算\(y\)轴上的宽度并，排序+扫一遍就行。
但是如何确定角度呢。
注意到，假如圆的相交关系不变，宽度的变化跟角度的变化是连续的。
先给出这个夹角的一个特解，然后找到每一段连续相交的圆的\(y\)最小和最大的两个，那么

\[w = \sum (rot(p_a, \theta + t).y - rot(p_b, \theta + t).y) + 2 \]

\[w = \sum (A.xsin\theta + A.ycos\theta - B.xsin\theta - B.ycos\theta) + 2 \]

求个导

\[dw = \sum (A.x cos\theta - A.ysin\theta - B.x cos \theta + B.y sin \theta) d\theta \]

\[dw = \sum (A.x - B.x) cos\theta - \sum (A.y - B.y)sin \theta) d\theta \]

令导数为\(0\)，得到\(tan\theta = \frac{\sum (A.y - B.y)}{\sum (A.x - B.x)}\)，于是可以得到\(\theta\)。

那么如何得到所有的特解呢，我们发现，两个圆投影是否相交当且仅当在两圆公切线处发生改变，于是我们枚举所有公切线角度，排序后就能得到每个区间了，可以取区间角平分线作为特解。

代码比较冗长。

#include <bits/stdc++.h>
#define pt(x) cout << x << endl;
#define Mid ((l + r) / 2)
#define lson (rt << 1)
#define rson (rt << 1 | 1)
using namespace std;
const double Pi = acos(-1.0);
const double eps = 1e-11;
struct Point {
	double x, y;
	Point() : x(0), y(0) {}
	Point(double x, double y) : x(x), y(y) {}
	Point operator+(const Point &a)const{ return Point(x + a.x, y + a.y); }
	Point operator-(const Point &a)const{ return Point(x - a.x, y - a.y); }
} a[109];
typedef Point Vector;
struct line{
	Point a, b;
	line() {}
	line(const Point &a,const Point &b) : a(a), b(b) {}
};
 
struct circle{
	Point c;
	double r;
	circle() : c(Point(0,0)), r(0) {}
	circle(const Point &c, double r) : c(c), r(r) {}
};
Point operator*(double c, const Point &a){ return Point(c * a.x, c * a.y); }
double abs(const Point &a){ return sqrt(a.x * a.x + a.y * a.y); }
Point rot(const Point &a, double theta){ return Point(a.x * cos(theta) - a.y * sin(theta), a.x * sin(theta) + a.y * cos(theta)); }
double arg(const Point &a){
	double t = atan2(a.y, a.x);
	return t < 0 ? t + 2 * Pi:t;
}
int tangent(const circle &C1,const circle &C2,vector<line> &res){
	double d = abs(C1.c - C2.c);
	if(d < eps) return 0;
 
	int c=0;
	if(C1.r + C2.r < d - eps){
		double t = acos((C1.r + C2.r) / d);
		res.push_back(line(C1.c + rot(C1.r / d * (C2.c - C1.c), t), C2.c + rot(C2.r / d * (C1.c - C2.c), t)));
		res.push_back(line(C1.c + rot(C1.r / d * (C2.c - C1.c),-t), C2.c + rot(C2.r / d * (C1.c - C2.c),-t)));
		c += 2;
	} else if(C1.r + C2.r < d + eps){
		Point p = C1.c + C1.r / d * (C2.c - C1.c);
		res.push_back(line(p, p + rot(C2.c - C1.c, Pi / 2)));
		c++;
	}
 
	if(abs(C1.r - C2.r) < d - eps){
		double t1 = acos((C1.r - C2.r) / d), t2 = Pi - t1;
		res.push_back(line(C1.c + rot(C1.r / d * (C2.c - C1.c), t1), C2.c + rot(C2.r / d * (C1.c - C2.c),-t2)));
		res.push_back(line(C1.c + rot(C1.r / d * (C2.c - C1.c),-t1), C2.c + rot(C2.r / d * (C1.c - C2.c), t2)));
		c+=2;
	} else if(abs(C1.r - C2.r) < d + eps){
		Point p = C1.c + C1.r / d * (C2.c - C1.c);
		res.push_back(line(p, p + rot(C2.c - C1.c, Pi / 2)));
		c++;
	}
 
	return c;
}
 
int n;
double maxn, minn, sima, simi;
double calc(double theta) {
	double tt[109] = {0};
	for(int k = 1; k <= n; k++) {
		tt[k] = rot(a[k], theta).y;
	}
	sort(tt + 1, tt + 1 + n);
	double mr = 0, sum = 0;
	sum = 2.0;
	mr = tt[1] + 1.0;
	for(int k = 1; k <= n; k++) {
		if(tt[k] + 1.0 < mr) continue;
		if(tt[k] - 1.0 < mr) {
			sum += tt[k] + 1.0 - mr;
		} else {
			sum += 2.0;
		}
		mr = tt[k] + 1.0;
	}
	return sum;
}
double calc_peak(double phi) {
	pair<double, int> p[109];
	for(int i = 1; i <= n; i++) {
		p[i] = make_pair(rot(a[i], phi).y, i);
	}
	sort(p + 1, p + 1 + n);
	int pre = 1;
	double x = 0, y = 0;
	for(int i = 2; i <= n + 1; i++) {
		if(i == n + 1 || p[i - 1].first + 1 < p[i].first - 1) {
			x += (a[p[i - 1].second].x - a[p[pre].second].x);
			y += (a[p[i - 1].second].y - a[p[pre].second].y);
			pre = i;
		}
	}
	double ans = arg({y, x});
	if(ans > Pi) ans -= Pi;
	return ans;
}
 
void work() {
	if(n == 0) exit(0);
	for(int i = 1; i <= n; i++) 
		scanf("%lf%lf", &a[i].x, &a[i].y);
	maxn = -1e18;
	minn = 1e18;
	vector<line> v;
	for(int i = 1; i <= n; i++) {
		for(int j = i + 1; j <= n; j++) {
			tangent(circle(a[i], 1.0), circle(a[j], 1.0), v);
		}
	}
	vector<double> th;
	for(int i = 0; i < v.size(); i++) {
		line &l = v[i];
		double tt = arg(l.b - l.a);
		if(tt < Pi) tt = Pi - tt;
		else tt = 2 * Pi - tt;
		th.push_back(tt);
	}
	th.push_back(0);
	th.push_back(Pi);
	sort(th.begin(), th.end());
	int tot = 0;
	for(int i = 0; i < th.size(); i++) 
		if(i == 0 || th[i] > th[i - 1] + eps) 
			th[tot++] = th[i];
	th.erase(th.begin() + tot, th.end());
	for(int i = 0; i < tot; i++) {
		double sum = calc(th[i]);
		if(sum > maxn + eps) {
			maxn = sum;
			sima = th[i];
		}
		if(sum < minn - eps) {
			minn = sum;
			simi = th[i];
		}
		if(i - 1 >= 0) {
			double phi0 = (th[i - 1] + th[i]) / 2;
			double phi1 = calc_peak(phi0);
			double sum = calc(phi1);
			if(sum > maxn + eps) {
				maxn = sum;
				sima = phi1;
			}
			if(sum < minn - eps) {
				minn = sum;
				simi = phi1;
			}
		}
	}
	printf("%.15f\n%.15f\n", simi, sima);
}
signed main()
{
	ios :: sync_with_stdio(0);
	cin.tie(0);
	while(~scanf("%d", &n)) work(); 
	return 0;
}