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

//H题
#include <complex>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <cstdlib>
#include <vector>
#include <algorithm>
#include <queue>
using namespace std;

const int N = 500;
const int TIM = 5000;
const double pi  = acos(-1.0);
const double eps = 1e-8;
typedef complex<double> Point;

double Det(const Point & a, const Point & b){
    return (conj(a) * b).imag();
}

double Dot(const Point & a, const Point & b){
    return (conj(a) * b).real();
}

int sgn(double x){
    if(abs(x) < eps) return 0;
    if(x < 0) return -1;
    return 1;
}

struct Line : public vector <Point>{
    Line(){};
    Line(const Point & a, const Point & b){
        push_back(a);
        push_back(b);
        return;
    }
};

Point Vec(const Line & a){
    return a[1] - a[0];
}

Point LineIntersection(const Line & a, const Line & b){
    Point u = a[0] - b[0];
    double k1 = Det(Vec(a), Vec(b));
    double k2 = Det(Vec(b), u);
    if(!sgn(k1))
        return a[0];
    return a[0] + Vec(a) * k2 / k1;
}

bool SegmentIntersection(const Line & a, const Line & b, int f){ 
    if(max(a[0].real(), a[1].real()) < min(b[0].real(), b[1].real())) return 0; 
    if(max(b[0].real(), b[1].real()) < min(a[0].real(), a[1].real())) return 0;
    if(max(a[0].imag(), a[1].imag()) < min(b[0].imag(), b[1].imag())) return 0;
    if(max(b[0].imag(), b[1].imag()) < min(a[0].imag(), a[1].imag())) return 0;
    Point t1 = Vec(a), t2 = Vec(b);
    if(f)
        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); 
    else
        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);  
}

int LSFLAG = 0;
Point LSIntersection(const Line & a, const Line & b){
    Point c = LineIntersection (a, b);
    LSFLAG = 0;
    if(sgn(Det(a[0] - b[0], Vec(b))) * sgn(Det(a[1] - b[0], Vec(b))) > 0 ){
        LSFLAG = -1;
        return c;
    }
    if(!sgn(Det(Vec(a), Vec(b)))) LSFLAG = -1;
    else{
        if(sgn(Dot(a[0] - c, a[1] - c)) <= 0) return c;
        else LSFLAG = -1;
    }
    return c;
}

Point CrossVP(const Line & v, const Point & u){
    if(!sgn(v[0].real() - v[1].real()) && !sgn(v[0].imag() - v[1].imag())) return v[0];
    //if(sgn(Dot(Vec(v), u - v[0])) < 0) return v[0];//删掉两行是直线 一行不删是线段 
    //if(sgn(Dot(Vec(v), u - v[1])) > 0) return v[1];//删掉这行是射线
    double a = Dot(Vec(v), u - v[0]);
    return v[0] + a / norm(Vec(v)) * Vec(v);
}

bool Cmp(const Point & a, const Point & b){
    if(! sgn(a.real() - b.real())){
        if(sgn(b.imag() - a.imag())) return 1;
        return 0;
    }
    if(sgn(b.real() - a.real()) > 0) return 1;
    return 0; 
}

Point con[N], p[N], ori;
int cn;

bool GetConvexHull(int n){
    //if(n < 3) return 0;
    sort(p + 1, p + n + 1, Cmp);
    cn = 0;
    con[++ cn] = p[1];
    con[++ cn] = p[2];
    for(int i = 3; i <= n; i ++){
        while(cn > 1 && sgn(Det(p[i] - con[cn - 1], con[cn] - con[cn - 1])) >= 0)
            cn --;
        con[++ cn] = p[i];
    }
    
    int num = cn;
    con[++ cn] = p[n - 1];
    for(int i = n -2; i; i --){
        while(cn > num && sgn(Det(p[i] - con[cn - 1], con[cn] - con[cn - 1])) >= 0) 
        cn --;
        con[++ cn] = p[i];
    }
    
    if(cn <= 3) return 0;
    return 1;
}

int T, Tcase, n;
double D;

Point LB(const Point & a, const Point & b){
    if(sgn(a.imag() - b.imag()) < 0) return a;
    else if(sgn(a.imag() - b.imag()) == 0 && sgn(a.real() - b.real()) <= 0) return a;
    return b;
}

Point RT(const Point & a, const Point & b){
    if(sgn(a.imag() - b.imag()) > 0) return a;
    else if(sgn(a.imag() - b.imag()) == 0 && sgn(a.real() - b.real()) >= 0) return a;
    return b;
}

bool Eq(const Point & a, const Point & b){
    if(!sgn(a.real() - b.real()) && !sgn(a.imag() - b.imag())) return 1;
    return 0;
}

void Solve(){
    scanf("%d%lf", &n, &D);
    for(int i = 1; i <= n; i ++){
        double x, y;
        scanf("%lf%lf", &x, &y);
        p[i] = (Point) {x, y};
    }
    
    double res = 0.0, d = pi / TIM;
    GetConvexHull(n);
    
    Point A, B, C;
    Line cur, cc;
    int cnt = 0;
    for(double th = 0; th < pi; th += d){
        cur = (Line) {(Point) {0, 0}, (Point){cos(th), sin(th)}};
        for(int i = 1; i <= cn; i ++){
            C = CrossVP(cur, con[i]);
            if(i == 1) A = C, B = C;
            else A = LB(A, C) , B = RT(B, C);
        }
        
        res += abs(A - B);
    }
    
    res = res / TIM / D;
        
    printf("Case #%d: %.4lf\n", ++Tcase, res);
    return;
}

int main(){
    scanf("%d", &T);
    while(T --)
        Solve();
    return 0;
}

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