计算几何学习4
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了计算几何学习4相关的知识，希望对你有一定的参考价值。
今天完成的内容很少
学习了一点半平面交
n^2的做法还是很平易近人
刚开始初始化一个大有界的平面 依次用每条直线去切割平面即可
原有的点如果在当前直线左侧一定会被保留
而原有多边形的线段 可能会在线段中间出现交点
在判断一下即可
不想加入重复的点 就在交点求出后判断一下
模板题 因为没注意题目的读入顺序可能顺时针可能逆时针
并且多边形也不一定凸 调了很久
其实 n^2的 HPI部分没错的
至于nlogn的照着板子敲了一发
WA了 刘汝佳老师板子很简洁
但是排序很迷 说是直线的极角排序 但没给重载
xiaotiantang学长的板子又很复杂
加上他的重载也没过
过程类似凸包 但是细节很多 还得好好再想想
至于今天的训练赛
失误就是今天的 A B 和 H
A 写了一发暴力 但是没加剪枝 纯二进制枚举 T了 就叫队友去写正解了
其实有人暴力过了
B 各种理解题意 开始以为是乘 后来又以为是只有 0 1 2
最后手算了一发 才知道是加
H 在比赛的时候有点石乐志
zkc告诉我可以枚举直线的夹角
再去暴力的算凸包在该斜率下的最远直线距离
但是在最后求距离的部分写的很繁琐
一直在想着用対踵点 但是没调好
赛后想了想 可以枚举垂线 这样距离就是凸包包上的点在垂线上的最远距离
然后就 sigma(dis) / D / 枚举次数即可
（还是因为一个小括号调了一阵
正解看这里
蒲丰抛针问题 长度为L的线段 随机投掷到无限组间距为D的平行直线中 与直线相交的概率为 2*L/(PI * D)
凸多边形为 C / (PI * D) 其中 C为周长
贴一发代码假装没有划水吧
H题
1 //H题 2 #include <complex> 3 #include <cstdio> 4 #include <cstring> 5 #include <cmath> 6 #include <iostream> 7 #include <cstdlib> 8 #include <vector> 9 #include <algorithm> 10 #include <queue> 11 using namespace std; 12 13 const int N = 500; 14 const int TIM = 5000; 15 const double pi = acos(-1.0); 16 const double eps = 1e-8; 17 typedef complex<double> Point; 18 19 double Det(const Point & a, const Point & b){ 20 return (conj(a) * b).imag(); 21 } 22 23 double Dot(const Point & a, const Point & b){ 24 return (conj(a) * b).real(); 25 } 26 27 int sgn(double x){ 28 if(abs(x) < eps) return 0; 29 if(x < 0) return -1; 30 return 1; 31 } 32 33 struct Line : public vector <Point>{ 34 Line(){}; 35 Line(const Point & a, const Point & b){ 36 push_back(a); 37 push_back(b); 38 return; 39 } 40 }; 41 42 Point Vec(const Line & a){ 43 return a[1] - a[0]; 44 } 45 46 Point LineIntersection(const Line & a, const Line & b){ 47 Point u = a[0] - b[0]; 48 double k1 = Det(Vec(a), Vec(b)); 49 double k2 = Det(Vec(b), u); 50 if(!sgn(k1)) 51 return a[0]; 52 return a[0] + Vec(a) * k2 / k1; 53 } 54 55 bool SegmentIntersection(const Line & a, const Line & b, int f){ 56 if(max(a[0].real(), a[1].real()) < min(b[0].real(), b[1].real())) return 0; 57 if(max(b[0].real(), b[1].real()) < min(a[0].real(), a[1].real())) return 0; 58 if(max(a[0].imag(), a[1].imag()) < min(b[0].imag(), b[1].imag())) return 0; 59 if(max(b[0].imag(), b[1].imag()) < min(a[0].imag(), a[1].imag())) return 0; 60 Point t1 = Vec(a), t2 = Vec(b); 61 if(f) 62 return (sgn(Det(a[0] - b[0], t2)) * sgn(Det(a[1] - b[0], t2)) <= 0) && (sgn(Det(b[0] - a[0], t1)) * sgn(Det(b[1] - a[0], t1)) <= 0); 63 else 64 return (sgn(Det(a[0] - b[0], t2)) * sgn(Det(a[1] - b[0], t2)) < 0) && (sgn(Det(b[0] - a[0], t1)) * sgn(Det(b[1] - a[0], t1)) < 0); 65 } 66 67 int LSFLAG = 0; 68 Point LSIntersection(const Line & a, const Line & b){ 69 Point c = LineIntersection (a, b); 70 LSFLAG = 0; 71 if(sgn(Det(a[0] - b[0], Vec(b))) * sgn(Det(a[1] - b[0], Vec(b))) > 0 ){ 72 LSFLAG = -1; 73 return c; 74 } 75 if(!sgn(Det(Vec(a), Vec(b)))) LSFLAG = -1; 76 else{ 77 if(sgn(Dot(a[0] - c, a[1] - c)) <= 0) return c; 78 else LSFLAG = -1; 79 } 80 return c; 81 } 82 83 Point CrossVP(const Line & v, const Point & u){ 84 if(!sgn(v[0].real() - v[1].real()) && !sgn(v[0].imag() - v[1].imag())) return v[0]; 85 //if(sgn(Dot(Vec(v), u - v[0])) < 0) return v[0];//删掉两行是直线 一行不删是线段 86 //if(sgn(Dot(Vec(v), u - v[1])) > 0) return v[1];//删掉这行是射线 87 double a = Dot(Vec(v), u - v[0]); 88 return v[0] + a / norm(Vec(v)) * Vec(v); 89 } 90 91 bool Cmp(const Point & a, const Point & b){ 92 if(! sgn(a.real() - b.real())){ 93 if(sgn(b.imag() - a.imag())) return 1; 94 return 0; 95 } 96 if(sgn(b.real() - a.real()) > 0) return 1; 97 return 0; 98 } 99 100 Point con[N], p[N], ori; 101 int cn; 102 103 bool GetConvexHull(int n){ 104 //if(n < 3) return 0; 105 sort(p + 1, p + n + 1, Cmp); 106 cn = 0; 107 con[++ cn] = p[1]; 108 con[++ cn] = p[2]; 109 for(int i = 3; i <= n; i ++){ 110 while(cn > 1 && sgn(Det(p[i] - con[cn - 1], con[cn] - con[cn - 1])) >= 0) 111 cn --; 112 con[++ cn] = p[i]; 113 } 114 115 int num = cn; 116 con[++ cn] = p[n - 1]; 117 for(int i = n -2; i; i --){ 118 while(cn > num && sgn(Det(p[i] - con[cn - 1], con[cn] - con[cn - 1])) >= 0) 119 cn --; 120 con[++ cn] = p[i]; 121 } 122 123 if(cn <= 3) return 0; 124 return 1; 125 } 126 127 int T, Tcase, n; 128 double D; 129 130 Point LB(const Point & a, const Point & b){ 131 if(sgn(a.imag() - b.imag()) < 0) return a; 132 else if(sgn(a.imag() - b.imag()) == 0 && sgn(a.real() - b.real()) <= 0) return a; 133 return b; 134 } 135 136 Point RT(const Point & a, const Point & b){ 137 if(sgn(a.imag() - b.imag()) > 0) return a; 138 else if(sgn(a.imag() - b.imag()) == 0 && sgn(a.real() - b.real()) >= 0) return a; 139 return b; 140 } 141 142 bool Eq(const Point & a, const Point & b){ 143 if(!sgn(a.real() - b.real()) && !sgn(a.imag() - b.imag())) return 1; 144 return 0; 145 } 146 147 void Solve(){ 148 scanf("%d%lf", &n, &D); 149 for(int i = 1; i <= n; i ++){ 150 double x, y; 151 scanf("%lf%lf", &x, &y); 152 p[i] = (Point) {x, y}; 153 } 154 155 double res = 0.0, d = pi / TIM; 156 GetConvexHull(n); 157 158 Point A, B, C; 159 Line cur, cc; 160 int cnt = 0; 161 for(double th = 0; th < pi; th += d){ 162 cur = (Line) {(Point) {0, 0}, (Point){cos(th), sin(th)}}; 163 for(int i = 1; i <= cn; i ++){ 164 C = CrossVP(cur, con[i]); 165 if(i == 1) A = C, B = C; 166 else A = LB(A, C) , B = RT(B, C); 167 } 168 169 res += abs(A - B); 170 } 171 172 res = res / TIM / D; 173 174 printf("Case #%d: %.4lf\n", ++Tcase, res); 175 return; 176 } 177 178 int main(){ 179 scanf("%d", &T); 180 while(T --) 181 Solve(); 182 return 0; 183 }
n^2的半平面交部分
int T, n, len; Point poly[N], pl[N]; Line l[N]; void CutPolygon(const Line & a, Point * poly, int &len){ if(!len) return; Point C, D, cc; int n = 0; for(int i = 1; i <= len; i ++){ int j = (i + 1); if(j > len) j -= len; C = poly[i]; D = poly[j]; if(sgn(Det(Vec(a), C - a[0])) >= 0) pl[++ n] = C; if(sgn(Det(Vec(a), C - D))){ cc = LSIntersection((Line) {C, D}, a); if(!LSFLAG && !Eq(C, cc) && !Eq(D, cc)) pl[++ n] = cc; } } len = n; for(int i = 1; i <= n; i ++) poly[i] = pl[i]; return; } void HPI(){ len = 0; poly[ ++ len] = (Point) {-inf, -inf}; poly[ ++ len] = (Point) {inf, -inf}; poly[ ++ len] = (Point) {inf, inf}; poly[ ++ len] = (Point) {-inf, inf}; for(int i = 1; i <= n; i ++){ int j = i + 1; if(j > n) j -= n; l[i] = (Line) {p[i], p[j]}; } for(int i = 1; i <= n; i ++) CutPolygon(l[i], poly, len); }
nlogn的明天再搞吧
明天是要搞好nlogn的半平面交 再做些相关习题吧
以上是关于计算几何学习4的主要内容，如果未能解决你的问题，请参考以下文章