码迷,mamicode.com
首页 > 其他好文 > 详细

计算几何学习4

时间:2017-07-21 00:03:10      阅读:252      评论:0      收藏:0      [点我收藏+]

标签:ase   顺时针   aot   name   cti   顺序   计算几何   main   目的   

今天完成的内容很少

学习了一点半平面交

 

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的半平面交 再做些相关习题吧

计算几何学习4

标签:ase   顺时针   aot   name   cti   顺序   计算几何   main   目的   

原文地址:http://www.cnblogs.com/drzdk/p/7214615.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!