博客
关于我
强烈建议你试试无所不能的chatGPT，快点击我
和项目组研究计算几何
阅读量：5972 次
发布时间：2019-06-19

本文共 15784 字，大约阅读时间需要 52 分钟。

一，凸包

二维凸包：可用极角排序后每次倍增一个点看是否在栈顶两个点内侧，在内侧就加入栈，否则弹出一个。优化是在省去了计算极角的计算量，采用xy排序，从x最小倍增一遍求上边界，再从x最大反过来求下边界。二维凸包构造想法比较简单。

凸多边形有一些有趣的性质

 

设边数为n则

内角和 =（n-2)×180°

 

外角和 = 360°

 

对角线的条数=C(n，2)-n=n(n-3)/2

欧拉公式 凸多边形有n个点，m条边，r个面，则 n - m + r = 2

 

 

1.POJ1113 注意精度，这里有个公式 城堡围墙长度最小值 = 城堡顶点坐标构成的散点集的凸包总边长 + 半径为L的圆周长

#include
     
      #include 
      
       #include
       
        using namespace std;#define eps 1e-8const double PI = acos(-1.0);struct point {    double x, y;    point() {    }    point(int a, int b) {        x = a, y = b;    }};int sgn(double a) {    if (fabs(a) < eps)        return 0;    if (a > eps)        return 1;    return -1;}double dis(point a, point b) {    return sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));}double cross(point op, point sp, point ep) {    return (sp.x - op.x) * (ep.y - op.y) - (ep.x - op.x) * (sp.y - op.y);}int cmp(point a, point b) {    return (a.y == b.y) ? (a.x < b.x) : a.y < b.y;}double graham(point pnt[], int n, point res[]) {    int i, len, top = 1;    sort(pnt, pnt + n, cmp);    res[0] = pnt[0];    res[1] = pnt[1];    for (i = 2; i < n; i++) {        while (top && sgn(cross(pnt[i], res[top], res[top - 1])) >= 0)            top--;        res[++top] = pnt[i];    }    len = top;    for (i = n - 2; i >= 0; i--) {        while (top != len && sgn(cross(pnt[i], res[top], res[top - 1])) >= 0)            top--;        res[++top] = pnt[i];    }    res[top] = res[0];    double ans = 0;    for (int p = 0; p < top; p++) {        ans += dis(res[p], res[p + 1]);    }    return ans;}int main() {//    freopen("data3.txt", "r", stdin);    point pnt[50003], res[50003];    int i, n;    double L;    while (scanf("%d%lf", &n, &L) != EOF) {        if (n == 2) {            for (i = 0; i < 2; i++)                scanf("%lf%lf", &pnt[i].x, &pnt[i].y);            printf("%.0lf\n", dis(pnt[0], pnt[1]) + L * 2 * PI);            continue;        }        for (i = 0; i < n; i++)            scanf("%lf%lf", &pnt[i].x, &pnt[i].y);        printf("%.0lf\n", graham(pnt, n, res) + L * 2 * PI);    }    return 0;}

 2.HDU 1872 枚举状态，使用位运算

#include 
     
      #include 
      
       #include 
       
        #include 
        
         #include 
         
          using namespace std;#define MAXN 20int N;#define eps 1e-8const double PI = acos(-1.0);struct point {    double x, y;    point() {    }    point(double a, double b) {        x = a, y = b;    }};struct Node {    point pos;    int val, len;} tree[MAXN];int sgn(double a) {    if (fabs(a) < eps)        return 0;    if (a > eps)        return 1;    return -1;}double dis(point a, point b) {    return sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));}double cross(point op, point sp, point ep) {    return (sp.x - op.x) * (ep.y - op.y) - (ep.x - op.x) * (sp.y - op.y);}int cmp(point a, point b) {    if (a.y < b.y)        return 1;    if (a.y == b.y)        if (a.x < b.x)            return 1;    return 0;}double graham(point pnt[], int n) {    point res[MAXN];    int i, len, top = 1;    sort(pnt, pnt + n, cmp);    res[0] = pnt[0];    res[1] = pnt[1];    for (i = 2; i < n; i++) {        while (top && sgn(cross(pnt[i], res[top], res[top - 1])) >= 0)            top--;        res[++top] = pnt[i];    }    len = top;    for (i = n - 2; i >= 0; i--) {        while (top != len && sgn(cross(pnt[i], res[top], res[top - 1])) >= 0)            top--;        res[++top] = pnt[i];    }    res[top] = res[0];    double ans = 0;    for (int p = 0; p < top; p++) {        ans += dis(res[p], res[p + 1]);    }    return ans;}struct Ans {    int bit;    double ex_len;    Ans() {        bit = 0, ex_len = 0;    }};int main() {//    freopen("data3.txt", "r", stdin);    point pnt[MAXN];    int tot;    int min_val, min_num;    int cas = 1;    while (scanf("%d", &N) && N) {        if (cas > 1)            puts("");        printf("Forest %d\n", cas++);        Ans ans;        min_val = min_num = 0x3fffffff;        for (int i = 0; i < N; i++) {            scanf("%lf%lf%d%d", &tree[i].pos.x, &tree[i].pos.y, &tree[i].val,                    &tree[i].len);        }        int cut_len, cut_val;        for (int bit = 0; bit < (1 << N); bit++) {            cut_len = cut_val = tot = 0;            for (int j = 0; j < N; j++) {                if (bit & (1 << j)) {                    cut_len += tree[j].len;                    cut_val += tree[j].val;                } else {                    pnt[tot].x = tree[j].pos.x, pnt[tot++].y = tree[j].pos.y;                }            }            if (cut_val > min_val)                continue;            double need_len = graham(pnt, tot);            if (need_len <= cut_len) {                if (min_val > cut_val                        || (min_val == cut_val && N - tot < min_num)) {                    min_val = cut_val, min_num = N - tot;                    ans.bit = bit;                    ans.ex_len = cut_len - need_len;                }            }        }        printf("Cut these trees:");        for (int i = 0; i < N; i++)            if (ans.bit & 1 << i)                printf(" %d", i + 1);        puts("");        printf("Extra wood: %.2lf\n", ans.ex_len);    }    return 0;}

 

三维凸包 倍增思想

#include
     
      #include
      
       #include
       
        #include
        
         using namespace std;const int MAXN = 505;const double eps = 1e-8;int zero(double x) {    if (fabs(x) < eps)        return 0;    return (x > eps) ? 1 : -1;}struct Point {    double x, y, z;    Point() {    }    Point(double xx, double yy, double zz) :            x(xx), y(yy), z(zz) {    }    Point operator -(const Point p1) {  //两向量之差        return Point(x - p1.x, y - p1.y, z - p1.z);    }    Point operator *(Point p) {   //叉乘        return Point(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x);    }    double operator ^(Point p) { //点乘        return (x * p.x + y * p.y + z * p.z);    }};struct CH3D {    struct face {        int a, b, c; //凸包一个面上三个点的编号        bool ok;  //表示该面是否属于最终凸包中的面    };    int n;   //顶点数    Point P[MAXN];  //初始顶点    int triangle_num;  //凸包表面的三角形数    face F[8 * MAXN];    int g[MAXN][MAXN]; //凸包表面的三角形    double vlen(Point a) {        return sqrt(a.x * a.x + a.y * a.y + a.z * a.z);    }    Point cross(const Point &a, const Point &b, const Point &c) {        return Point((b.y - a.y) * (c.z - a.z) - (b.z - a.z) * (c.y - a.y),                -((b.x - a.x) * (c.z - a.z) - (b.z - a.z) * (c.x - a.x)),                (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x));    }    //三角形面积*2    double area(Point a, Point b, Point c) {        return vlen((b - a) * (c - a));    }    //四面体a到其他三点有向距离有向体积*6    double volume(Point a, Point b, Point c, Point d) {        return (b - a) * (c - a) ^ (d - a);    }    //点在面同向返回正数    int dblcmp(Point &p, face &f) {        double v = volume(P[f.a], P[f.b], P[f.c], p);        if (fabs(v) < eps)            return 0;        return (v > eps) ? 1 : -1;    }    void deal(int p, int a, int b) {        int f = g[a][b];        face add;        if (F[f].ok) {            if (dblcmp(P[p], F[f]) > 0)                dfs(p, f);            else {                add.a = b;                add.b = a;                add.c = p;                add.ok = 1;                g[p][b] = g[a][p] = g[b][a] = triangle_num;                F[triangle_num++] = add;            }        }    }    void dfs(int p, int now) {        F[now].ok = 0;        deal(p, F[now].b, F[now].a);        deal(p, F[now].c, F[now].b);        deal(p, F[now].a, F[now].c);    }    //调整使前四个点不共面    bool check_4point() {        int i, flag = 1;        for (i = 1; i < n; i++) {            if (vlen(P[0] - P[i]) > eps) {                swap(P[1], P[i]);                flag = 0;                break;            }        }        if (flag)            return 0;        flag = 1;        for (i = 2; i < n; i++) {            if (vlen((P[0] - P[1]) * (P[1] - P[i])) > eps) {                swap(P[2], P[i]);                flag = 0;                break;            }        }        if (flag)            return 0;        flag = 1;        for (i = 3; i < n; i++) {            if (fabs((P[0] - P[1]) * (P[1] - P[2]) ^ (P[0] - P[i])) > eps) {                swap(P[3], P[i]);                return 1;            }        }        return 0;    }    //构建三维凸包    void build() {        int i, j, tmp;        face add;        triangle_num = 0;        if (n < 4)            return;        if (!check_4point()) //调整前4个点，找不到不共面4个点则无法构建凸包; 若以保证,则可去掉            return;        for (i = 0; i < 4; i++) { //根据前4个点产生4个三角面            add.a = (i + 1) % 4;            add.b = (i + 2) % 4;            add.c = (i + 3) % 4;            add.ok = 1;            if (dblcmp(P[i], add) > 0)                swap(add.b, add.c);            g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = triangle_num;            F[triangle_num++] = add;        }        for (i = 4; i < n; i++) {            for (j = 0; j < triangle_num; j++) {                if (F[j].ok && dblcmp(P[i], F[j]) > 0) {                    dfs(i, j);                    break;                }            }        }        tmp = triangle_num;        for (i = triangle_num = 0; i < tmp; i++)            if (F[i].ok) {                F[triangle_num++] = F[i];            }    }    //三维凸包表面积    double area() {        double res = 0.0;        if (n == 3) {            Point p = cross(P[0], P[1], P[2]);            res = vlen(p) / 2.0;            return res;        }        for (int i = 0; i < triangle_num; i++)            res += fabs(area(P[F[i].a], P[F[i].b], P[F[i].c]));        return res / 2.0;    }    //三维凸包体积    double volume() {        double res = 0.0;        Point tmp(0, 0, 0);        for (int i = 0; i < triangle_num; i++)            res += volume(tmp, P[F[i].a], P[F[i].b], P[F[i].c]);        return fabs(res / 6.0);    }    //判断两个三角形共面    bool same(int s, int t) {        Point a = P[F[s].a];        Point b = P[F[s].b];        Point c = P[F[s].c];        return !zero(volume(a, b, c, P[F[t].a]))                && !zero(volume(a, b, c, P[F[t].b]))                && !zero(volume(a, b, c, P[F[t].c]));    }    //三维凸包表面多边形个数    int polygon() {        int i, j, res, flag;        for (i = res = 0; i < triangle_num; i++) {            flag = 1;            for (j = 0; j < i; j++)                if (same(i, j)) {                    flag = 0;                    break;                }            res += flag;        }        return res;    }};CH3D hull;int main() {//    freopen("data2.txt", "r", stdin);    while (~scanf("%d", &hull.n)) {        for (int i = 0; i < hull.n; i++)            scanf("%lf%lf%lf", &hull.P[i].x, &hull.P[i].y, &hull.P[i].z);        hull.build();        printf("%d\n", hull.polygon());//        res = hull.area();//        printf("%.3lf\n", res);    }    return 0;}

 

 

 

二，半平面

对所有线极角排序去重，之后维护一个直线双向栈，倍增一条线时，检查双向栈两端栈顶两条线的交点是否在新线内侧，内侧就加入栈，否则弹出一条线，最后双向栈中的各条线组成图形的核，可以在这个区域里看到原图形的每一个位置。用这个求内切的最大的圆，就是把各边按二分半径收缩后看是否还存在核的最大半径即解。

 

#include 
     
      #include 
      
       #include 
       
        using namespace std;const double eps = 1e-8;const int maxn = 105;int deq[maxn], tail, head, order[maxn], ln ,pn;struct Point {    double x, y;} p[maxn];struct Line {    Point a, b;    double angle;} l[maxn];int sgn(double a){    if(fabs(a)
        
         eps) return 1;    return -1;}int cross(Point op, Point sp, Point ep){    return (sp.x-op.x) * (ep.y-op.y) - (ep.x-op.x) * (sp.y-op.y);}bool cmp(int u, int v) {    int d = sgn(l[u].angle-l[v].angle);    if (!d) return sgn(cross(l[u].a, l[v].a, l[v].b)) < 0; //极角相同时，只保留最靠里面的那条,将更靠半平面里面的放在前面去重时删除    return d < 0;}void getIntersect(Line l1, Line l2, Point& p) { //求两线交点    double dot1,dot2;    dot1 = cross(l2.a, l1.b, l1.a);    dot2 = cross(l1.b, l2.b, l1.a);    p.x = (l2.a.x * dot2 + l2.b.x * dot1) / (dot2 + dot1);    p.y = (l2.a.y * dot2 + l2.b.y * dot1) / (dot2 + dot1);}bool judge(Line l0, Line l1, Line l2) {    Point p;    getIntersect(l1, l2, p);    return sgn(cross(p, l0.a, l0.b)) > 0; //大于小于符号与上面cmp()中注释处相反}bool halfPlaneIntersection() {    int i, j;    for (i = 0; i < ln; i++) order[i] = i;    sort(order, order+ln, cmp);    for (i = 1, j = 0; i < ln; i++)        if (sgn(l[order[i]].angle-l[order[j]].angle) > 0)            order[++j] = order[i];    ln = j + 1;    deq[0] = order[0];    deq[1] = order[1];    head = 0;    tail = 1;    for (i = 2; i < ln; i++) {        while (head < tail && judge(l[order[i]], l[deq[tail-1]], l[deq[tail]])) tail--;        while (head < tail && judge(l[order[i]], l[deq[head+1]], l[deq[head]])) head++;        deq[++tail] = order[i];    }    while (head < tail && judge(l[deq[head]], l[deq[tail-1]], l[deq[tail]])) tail--;    while (head < tail && judge(l[deq[tail]], l[deq[head+1]], l[deq[head]])) head++;    if (head + 1 >= tail) return 0;    //求出核多边形点p[]    deq[++tail] = deq[head];    for(pn=0,i=head;i
         
          1) printf("\n");        printf("Floor #%d\n",cnt++);        for(i=0;i

 

 

三，三维几何 线段

与二维相似综合使用点乘叉乘计算，三维的叉乘比较特殊，两个向量的叉乘为垂直这两个向量的向量，根据这个性质适合求法线，四面体的高等等。和二维一样，点乘判垂直，叉乘判平行，点乘另一个应用是求一个向量在另一个单位向量上的投影长度。

HDU 4741

//三维几何函数库#include 
     
      #define eps 1e-8#define zero(x) (((x)>0?(x):-(x))
      
       eps&&        vlen(xmult(subt(p,s.b),subt(p,s.c)))>eps&&vlen(xmult(subt(p,s.c),subt(p,s.a)))>eps;}int dot_inplane_ex(point3 p,point3 s1,point3 s2,point3 s3){    return dot_inplane_in(p,s1,s2,s3)&&vlen(xmult(subt(p,s1),subt(p,s2)))>eps&&        vlen(xmult(subt(p,s2),subt(p,s3)))>eps&&vlen(xmult(subt(p,s3),subt(p,s1)))>eps;}//判两点在线段同侧,点在线段上返回0,不共面无意义int same_side(point3 p1,point3 p2,line3 l){    return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))>eps;}int same_side(point3 p1,point3 p2,point3 l1,point3 l2){    return dmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))>eps;}//判两点在线段异侧,点在线段上返回0,不共面无意义int opposite_side(point3 p1,point3 p2,line3 l){    return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))<-eps;}int opposite_side(point3 p1,point3 p2,point3 l1,point3 l2){    return dmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))<-eps;}//判两点在平面同侧,点在平面上返回0int same_side(point3 p1,point3 p2,plane3 s){    return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))>eps;}int same_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3){    return dmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))>eps;}//判两点在平面异侧,点在平面上返回0int opposite_side(point3 p1,point3 p2,plane3 s){    return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))<-eps;}int opposite_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3){    return dmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))<-eps;}//判两直线平行int parallel(line3 u,line3 v){    return vlen(xmult(subt(u.a,u.b),subt(v.a,v.b)))
       
        > T;    while (T--) {        point3 a[2], b[2];        scanf("%lf%lf%lf", &a[0].x, &a[0].y, &a[0].z);        scanf("%lf%lf%lf", &a[1].x, &a[1].y, &a[1].z);        scanf("%lf%lf%lf", &b[0].x, &b[0].y, &b[0].z);        scanf("%lf%lf%lf", &b[1].x, &b[1].y, &b[1].z);        point3 dir = xmult(subt(a[0], a[1]), subt(b[0], b[1])); //叉积与两直线垂直的向量        double dis = 0;        if (vlen(dir) > eps) { //判断共面            double d = 1 / vlen(dir);            dir.x = dir.x * d;            dir.y = dir.y * d;            dir.z = dir.z * d;            dis = dmult(subt(a[0], b[0]), dir); //与单位向量的点积为这个方向的长度            dir.x = dir.x * dis;            dir.y = dir.y * dis;            dir.z = dir.z * dis;        }        printf("%.6f\n", fabs(dis));        point3 pa = intersection(a[0], a[1], plusv(b[0], dir),plusv(b[1], dir));        point3 pb = subt(pa, dir);        printf("%lf %lf %lf %lf %lf %lf\n", pa.x, pa.y, pa.z, pb.x, pb.y, pb.z);    }    return 0;}

 

转载于:https://www.cnblogs.com/updateofsimon/p/3423730.html

你可能感兴趣的文章
linux的日志服务器关于屏蔽一些关键字的方法
查看>>
mysql多实例实例化数据库
查看>>
javascript 操作DOM元素样式
查看>>
HBase 笔记3
查看>>
【Linux】Linux 在线安装yum
查看>>
Atom 编辑器系列视频课程
查看>>
[原][osgearth]osgearthviewer读取earth文件，代码解析（earth文件读取的一帧）
查看>>
使用dotenv管理环境变量
查看>>
Vuex学习
查看>>
bootstrap - navbar
查看>>
FastDFS存储服务器部署
查看>>
Android — 创建和修改 Fragment 的方法及相关注意事项
查看>>
swift基础之_swift调用OC/OC调用swift
查看>>
ElasticSearch Client详解
查看>>
mybatis update返回值的意义
查看>>
expdp 详解及实例
查看>>
通过IP判断登录地址
查看>>
深入浅出JavaScript （五） 详解Document.write()方法
查看>>
Beta冲刺——day6
查看>>
在一个程序中调用另一个程序并且传输数据到选择屏幕执行这个程序
查看>>
喝酒易醉，品茶养心，人生如梦，品茶悟道，何以解忧？唯有杜康！-- 愿君每日到此一游！
当前时间: 2024-12-29 08:55:18   当前IP: 18.222.161.119   联系邮箱:javaeecc@qq.com     Copyright © 2020 - 2022 baihongyu.com 京ICP备2021015314号-2
强烈建议你试试无所不能的CHAT-GPT，快点击我
强烈建议你试试无所不能的CHAT-GPT，快点击我