zoukankan      html  css  js  c++  java
  • 计算几何学习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的半平面交 再做些相关习题吧

  • 相关阅读:
    前沿技术解密——VirtualDOM
    Ques核心思想——CSS Namespace
    Unix Pipes to Javascript Pipes
    Road to the future——伪MVVM库Q.js
    聊聊CSS postproccessors
    【译】十款性能最佳的压缩算法
    Kafka Streams开发入门(9)
    Kafka Streams开发入门(8)
    【译】Kafka Producer Sticky Partitioner
    【译】99th Percentile Latency at Scale with Apache Kafka
  • 原文地址:https://www.cnblogs.com/drzdk/p/7214615.html
Copyright © 2011-2022 走看看