hihoCoder 1231 Border Length(ACM-ICPC国际大学生程序设计竞赛北京赛区(2015)网络赛 E)

题目链接：
http://hihocoder.com/problemset/problem/1231
题意：依次给出一个简单多边形上的顶点坐标（鸟飞行的轨迹）和一个圆（未名湖），求多边形和圆的重合部分的周长。

之前做过一个求圆和简单多边形重合部分面积的题，用类似于求任意多边形面积的方法就可以搞定。

然而这个题的方法是求多边形所有边在圆里面的部分的长度+圆在多边形内的弧的长度。

对多边形，只需要求出多边形和圆的所有交点，用这些交点把多边形的每条边分成几部分，把在圆内部分的长度加上。

对于圆，还要用上面求出的交点把圆分成很多段弧，依次判断这些弧是不是在多边形内。（一段弧不是完全在多边形里面就是完全在多边形外面，所以只需要判断弧的中点是不是在多边形里面）对于交点个数<=1的情况（圆和多边形的某一条边相切于一点），多边形与圆的重合部分要么只有这个点，要么是整个圆。只需要判断圆心是否在多边形内，就可以判断结果是否包括了整个圆的周长。

关于判断点是否在简单多边形内，找到了两种方法：

1.用射线和多边形的交点判断。从要判断的点P任意作一条射线，如果这条射线和多边形的交点个数为奇数，则点P在多边形内，否则点P在多边形外。在这种算法中，对于点在多边形上（点P与多边形某个顶点重合或在多边形的某一条边上）应该特殊处理。（这段代码是抄的某个板，其实没看太懂= =）

2.回转数法 在平面内，封闭曲线绕一个固定点的回转数定义为曲线沿逆时针方向绕过该点的次数，是一个整数。如果曲线沿顺时针方向绕过该点则回转数为负数。显然，对于多边形外的一点，曲线绕该点的回转数为0。(https://en.wikipedia.org/wiki/Winding_number)

（一开始想用点到边的两个端点的连线的夹角*叉积的符号求，发现点在多边形上的时候跪了……其实沿某个固定方向求夹角之和就可以）

参考了一下这篇文章里的js代码：
http://www.tuicool.com/articles/RnUBf2
//写计算几何没有板真是坑爹啊……

贴代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <vector>
#include <assert.h>
using namespace std;

const double eps = 1e-8;
const double PI = acos(-1.0);

const int ONEDGE = 2;
const int INSIDE = 1;
const int OUTSIDE = 0;

int sgn(double a){
if (fabs(a) < eps) return 0;
return a < 0 ? -1 : 1;
}

class Point {
public:
double x, y;
Point(double _x = 0, double _y = 0) :x(_x), y(_y){}

double dist_to(Point a){
return sqrt((x - a.x)*(x - a.x) + (y - a.y)*(y - a.y));
}
};

typedef Point Vector;
Vector operator + (const Vector &a, const Vector &b){
return Vector(a.x + b.x, a.y + b.y);
}
Vector operator - (const Vector &a, const Vector &b){
return Vector(a.x - b.x, a.y - b.y);
}
Vector operator * (const Vector &a, const double k){
return Vector(a.x*k, a.y*k);
}
Vector operator * (const double k, const Vector &a){
return Vector(a.x*k, a.y*k);
}
double dmul(Vector a, Vector b){
return a.x*b.x + a.y*b.y;
}
double xmul(Vector a, Vector b){
return a.x*b.y - a.y*b.x;
}

class Circle{
public:
Point o; double r;
};
int is_point_in_circle(Circle c, Point p){
int flag = sgn(c.o.dist_to(p) - c.r);
if (flag > 0)return OUTSIDE;
if (flag == 0)return ONEDGE;
if (flag < 0)return INSIDE;
}
//求弧的中点
//从ang1逆时针旋转到arg2的弧，范围为[-PI, PI]
Point get_arc_mid(Circle c, double ang1, double ang2){
if (ang2 < ang1){
ang2 += 2 * PI;
}
double ang = (ang1 + ang2) / 2;
return c.o + c.r*Vector(cos(ang), sin(ang));
}
//求弧的长度
//从ang1逆时针旋转到arg2的弧，范围为[-PI, PI]
double get_arc_len(Circle c, double ang1, double ang2){
if (ang2 < ang1){
ang2 += 2 * PI;
}
double ang = fabs(ang1 - ang2);
return c.r*ang;
}

class Segment{
public:
Point a, b;
Segment(Point _a, Point _b) :a(_a), b(_b){};
//线段长度
double length() {
return a.dist_to(b);
}
//线段中点
Point mid() {
return Point((a.x + b.x) / 2, (a.y + b.y) / 2);
}
};
bool is_point_on_segment(Segment se, Point p){
bool ans = false;
if (sgn(p.dist_to(se.a)) == 0 || sgn(p.dist_to(se.b)) == 0)
ans = true;
else{
double minx = min(se.a.x, se.b.x), maxx = max(se.a.x, se.b.x);
double miny = min(se.a.y, se.b.y), maxy = max(se.a.y, se.b.y);
if (sgn(minx - p.x) <= 0 && sgn(p.x - maxx) <= 0 && sgn(miny - p.y) <= 0 && sgn(p.y - maxy) <= 0) {
/*
在这里被坑了
要判断线段到点的距离而不是直接判断三角形面积
有可能会出现线段特别长的情况，判断面积会出错
double xx = xmul(p - se.a, p - se.b);
double dis = se.length();
assert(sgn(xx - (xx / dis)) == 0);
*/
ans = (fabs(xmul(p - se.a, p - se.b) / se.length()) < eps);
}
}
return ans;
}
bool is_segment_in_circle(Circle ci, Segment se){
Point mid = se.mid();
return sgn(ci.r - ci.o.dist_to(mid)) > 0;
}

vector<Point> get_circle_segment_intersection(Circle ci, Segment se){

vector<Point> res;
Vector v = se.b - se.a;
Point p0 = se.a;
double _a = v.x, _b = p0.x - ci.o.x, _c = v.y, _d = p0.y - ci.o.y;
double a = _a*_a + _c*_c;
double b = 2 * (_a * _b + _c * _d);
double c = _b * _b + _d * _d - ci.r * ci.r;

double delta = b*b - 4 * a*c;
if(sgn(delta) > 0) {
double k1 = (-b + sqrt(delta)) / (2 * a);
double k2 = (-b - sqrt(delta)) / (2 * a);
Point p1 = se.a + k1*v, p2 = se.a + k2*v;
if (is_point_on_segment(se, p1))
res.push_back(p1);
if (is_point_on_segment(se, p2))
res.push_back(p2);
}
else if(sgn(delta) >= 0) {
double k1 = (-b) / (2 * a);
Point p1 = se.a + k1*v;
if (is_point_on_segment(se, p1))
res.push_back(p1);
}

return res;
}

class Polygon{
public:
vector<Point> p;
};
int is_point_in_polygon(Polygon polygon, Point a){
//三种方法都可以过
Point &p = a;

//winding number
double sum = 0.0, tmp;
for (int i = 0; i < polygon.p.size(); i++){
Point a = polygon.p[i], b = polygon.p[(i + 1) % polygon.p.size()];
Segment se(a, b);
if (is_point_on_segment(se, p)){
return true;
}
Vector v1 = a - p, v2 = b - p;
double angle = (atan2(v1.y, v1.x) - atan2(v2.y, v2.x));
if (angle >= (PI)){
angle = angle - (PI * 2);
}
else if (angle < (-PI)){
angle = angle + (PI * 2);
}
sum = sum + angle;
}
sum = fabs(sum);
//除PI或者(2*PI)貌似都可以
int i = (int)(sum / (2 * PI) + 0.5);
return i != 0;

/*
//ray casting
int i, j; bool ans = false;
for (i = 0, j = polygon.p.size() - 1; i < polygon.p.size(); j = i++){
if (((polygon.p[i].y>p.y) != (polygon.p[j].y > p.y)) &&
(p.x < (polygon.p[j].x - polygon.p[i].x) * (p.y - polygon.p[i].y) / (polygon.p[j].y - polygon.p[i].y) + polygon.p[i].x))
ans = !ans;
}
return ans;
*/
/*
//ray casting
Point q;
int i = 0;
int counter;
int len = polygon.p.size();
const double offset = 1e8;
while (i<len)
{
q.x = rand() + offset;//随机取一个足够远的点q
q.y = rand() + offset;//以p为起点q为终点做射线L
for (counter = i = 0; i<len; i++) //依次对多边形的每条边进行考察
{
Point t1 = polygon.p[i] - a;
Point t2 = polygon.p[(i + 1) % len] - a;
Point t3 = q - a;
Point t4 = polygon.p[i] - a;
Point t5 = a - polygon.p[i];
Point t6 = polygon.p[(i + 1) % len] - polygon.p[i];
Point t7 = q - polygon.p[i];
if (fabs(xmul(t1, t2))<eps &&
(polygon.p[i].x - a.x)*(polygon.p[(i + 1) % len].x - a.x)<eps && (polygon.p[i].y - a.y)*(polygon.p[(i + 1) % len].y - a.y)<eps)
return ONEDGE; //点p在边上,返回on_edge
else if (fabs(xmul(t3, t4))<eps) break; //点polygon.p[i]在射线pq上，停止本循环，另取q
else if (xmul(t4, t3)*xmul(t2, t3)<-eps &&
xmul(t5, t6)*xmul(t7, t6)<-eps)
counter++;
}
}
return (counter & 1) ? INSIDE : OUTSIDE;
*/
}

int n;
double ans = 0;
vector<double> div_angle;
Polygon p; Circle c;

//多边形在圆内的长度
void func1(){
vector<Segment> seg;
for (int i = 0; i < n; i++){
int j = (i + 1) % n;
seg.clear();
Point a = p.p[i], b = p.p[j];
Segment se(a, b);
vector<Point> inter_point = get_circle_segment_intersection(c, se);

if (inter_point.size() == 0){
seg.push_back(se);
}
else if (inter_point.size() == 1){
Point p1 = inter_point[0];
Segment s1 = Segment(a, p1), s2 = Segment(p1, b);
seg.push_back(s1); seg.push_back(s2);
}
else {
Point p1 = inter_point[0], p2 = inter_point[1];
if (a.dist_to(p1) > a.dist_to(p2)){
swap(p1, p2);
}
Segment s1 = Segment(a, p1), s2 = Segment(p1, p2), s3 = Segment(p2, b);
seg.push_back(s1); seg.push_back(s2); seg.push_back(s3);
}

for (int k = 0; k < seg.size(); k++){
if (is_segment_in_circle(c, seg[k]))
ans += seg[k].length();
}
for (int k = 0; k < inter_point.size(); k++){
Vector v = inter_point[k] - c.o;
div_angle.push_back(atan2(v.y, v.x));
}
}
}

//圆在多边形内的长度
void func2(){
sort(div_angle.begin(), div_angle.end());
div_angle.erase(unique(div_angle.begin(), div_angle.end()), div_angle.end());
if (div_angle.size() <= 1){
if (is_point_in_polygon(p, c.o)){
ans += 2 * PI*c.r;
}
}
else {
div_angle.push_back(div_angle[0] + 2 * PI);
int size = div_angle.size() - 1;
for (int i = 0; i < size; i++){
double a1 = div_angle[i], a2 = div_angle[i + 1];
Point mid = get_arc_mid(c, a1, a2);
if (is_point_in_polygon(p, mid)){
ans += get_arc_len(c, a1, a2);
}
}
}
}

int main(){

while (scanf("%d", &n) && n){
ans = 0;
div_angle.clear();
p.p.clear();

Point tmp;
for (int i = 0; i < n; i++){
scanf("%lf%lf", &tmp.x, &tmp.y);
p.p.push_back(tmp);
}
scanf("%lf%lf%lf", &c.o.x, &c.o.y, &c.r);

func1();
func2();

printf("%.0lf\n", ans);
}
}