计算几何学习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 }

View Code

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);
}

View Code

nlogn的明天再搞吧

明天是要搞好nlogn的半平面交再做些相关习题吧

posted on 2017-07-20 21:34 drzdk 阅读(133) 评论(0) 编辑收藏举报