POJ 2225 / ZOJ 1438 / UVA 1438 Asteroids --三维凸包,求多面体重心

题意: 两个凸多面体,可以任意摆放,最多贴着,问他们重心的最短距离。

解法: 由于给出的是凸多面体,先构出两个三维凸包,再求其重心,求重心仿照求三角形重心的方式,然后再求两个多面体的重心到每个多面体的各个面的最短距离,然后最短距离相加即为答案,因为显然贴着最优。

求三角形重心见此: http://www.cnblogs.com/whatbeg/p/4234518.html

代码:(模板借鉴网上模板)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <string>
#include <vector>
#include <set>
#define Mod 1000000007
#define eps 1e-8
#define lll __int64
#define ll long long
using namespace std;
#define N 100007
#define MAXV 505

//三维点
struct pt{
    double x, y, z;
    pt(){}
    pt(double _x, double _y, double _z): x(_x), y(_y), z(_z){}
    pt operator - (const pt p1){return pt(x - p1.x, y - p1.y, z - p1.z);}
    pt operator * (pt p){return pt(y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x);}        //叉乘
    double operator ^ (pt p){return x*p.x+y*p.y+z*p.z;}                            //点乘
};

//pt operator - (const pt p,const pt p1){return pt(p.x - p1.x, p.y - p1.y, p.z - p1.z);}
//pt operator ** (pt p,pt p1){return pt(p.y*p1.z-p.z*p1.y, p.z*p1.x-p.x*p1.z, p.x*p1.y-p.y*p1.x);}        //叉乘
//double operator ^^ (pt p1,pt p){return p1.x*p.x+p1.y*p.y+p1.z*p.z;}

struct _3DCH{
    struct fac{
        int a, b, c;    //表示凸包一个面上三个点的编号
        bool ok;        //表示该面是否属于最终凸包中的面
    };

    int n;    //初始点数
    pt P[MAXV];    //初始点

    int cnt;    //凸包表面的三角形数
    fac F[MAXV*8]; //凸包表面的三角形

    int to[MAXV][MAXV];

    double vlen(pt a){return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);}    //向量长度
    double area(pt a, pt b, pt c){return vlen((b-a)*(c-a));}    //三角形面积*2
    double volume(pt a, pt b, pt c, pt d){return (b-a)*(c-a)^(d-a);}    //四面体有向体积*6

    //正:点在面同向
    double ptof(pt &p, fac &f){
        pt m = P[f.b]-P[f.a], n = P[f.c]-P[f.a], t = p-P[f.a];
        return (m * n) ^ t;
    }
    pt pvec(fac s) {
        pt k1 = (P[s.a]-P[s.b]), k2 = (P[s.b]-P[s.c]);
        return (k1*k2);
    }
    double ptoplane(pt p,fac s){
        return fabs(pvec(s)^(p-P[s.a]))/vlen(pvec(s));
    }

    void deal(int p, int a, int b){
        int f = to[a][b];
        fac add;
        if (F[f].ok){
            if (ptof(P[p], F[f]) > eps)
                dfs(p, f);
            else{
                add.a = b, add.b = a, add.c = p, add.ok = 1;
                to[p][b] = to[a][p] = to[b][a] = cnt;
                F[cnt++] = add;
            }
        }
    }

    void dfs(int p, int cur){
        F[cur].ok = 0;
        deal(p, F[cur].b, F[cur].a);
        deal(p, F[cur].c, F[cur].b);
        deal(p, F[cur].a, F[cur].c);
    }

    bool same(int s, int t){
        pt &a = P[F[s].a], &b = P[F[s].b], &c = P[F[s].c];
        return fabs(volume(a, b, c, P[F[t].a])) < eps && fabs(volume(a, b, c, P[F[t].b])) < eps && fabs(volume(a, b, c, P[F[t].c])) < eps;
    }

    //构建三维凸包
    void construct(){
        cnt = 0;
        if (n < 4)
            return;

        /*********此段是为了保证前四个点不公面,若已保证,可去掉********/
        bool sb = 1;
        //使前两点不公点
        for (int i = 1; i < n; i++){
            if (vlen(P[0] - P[i]) > eps){
                swap(P[1], P[i]);
                sb = 0;
                break;
            }
        }
        if (sb)return;

        sb = 1;
        //使前三点不公线
        for (int i = 2; i < n; i++){
            if (vlen((P[0] - P[1]) * (P[1] - P[i])) > eps){
                swap(P[2], P[i]);
                sb = 0;
                break;
            }
        }
        if (sb)return;

        sb = 1;
        //使前四点不共面
        for (int 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]);
                sb = 0;
                break;
            }
        }
        if (sb)return;
        /*********此段是为了保证前四个点不公面********/

        fac add;
        for (int i = 0; i < 4; i++){
            add.a = (i+1)%4, add.b = (i+2)%4, add.c = (i+3)%4, add.ok = 1;
            if (ptof(P[i], add) > 0)
                swap(add.b, add.c);
            to[add.a][add.b] = to[add.b][add.c] = to[add.c][add.a] = cnt;
            F[cnt++] = add;
        }

        for (int i = 4; i < n; i++){
            for (int j = 0; j < cnt; j++){
                if (F[j].ok && ptof(P[i], F[j]) > eps){
                    dfs(i, j);
                    break;
                }
            }
        }
        int tmp = cnt;
        cnt = 0;
        for (int i = 0; i < tmp; i++){
            if (F[i].ok){
                F[cnt++] = F[i];
            }
        }
    }

    //表面积
    double area(){
        double ret = 0.0;
        for (int i = 0; i < cnt; i++){
            ret += area(P[F[i].a], P[F[i].b], P[F[i].c]);
        }
        return ret / 2.0;
    }

    //体积
    double volume(){
        pt O(0, 0, 0);
        double ret = 0.0;
        for (int i = 0; i < cnt; i++) {
            ret += volume(O, P[F[i].a], P[F[i].b], P[F[i].c]);
        }
        return fabs(ret / 6.0);
    }

    pt BaryCenter() {
        pt O(0, 0, 0);
        double ret = 0.0,sumvolume = 0.0, sumx = 0.0, sumy = 0.0, sumz = 0.0;
        for(int i=0;i<cnt;i++) {
            double Vol = volume(O, P[F[i].a], P[F[i].b], P[F[i].c]);
            sumvolume += Vol;
            sumx += (P[F[i].a].x + P[F[i].b].x + P[F[i].c].x)*Vol;
            sumy += (P[F[i].a].y + P[F[i].b].y + P[F[i].c].y)*Vol;
            sumz += (P[F[i].a].z + P[F[i].b].z + P[F[i].c].z)*Vol;
        }
        return pt(sumx/sumvolume/4, sumy/sumvolume/4, sumz/sumvolume/4);
    }

    //表面三角形数
    int facetCnt_tri(){
        return cnt;
    }

    //表面多边形数
    int facetCnt(){
        int ans = 0;
        for (int i = 0; i < cnt; i++){
            bool nb = 1;
            for (int j = 0; j < i; j++){
                if (same(i, j)){
                    nb = 0;
                    break;
                }
            }
            ans += nb;
        }
        return ans;
    }
};

_3DCH hull,hull2;    //内有大数组,不易放在函数内

int main()
{
    while (scanf("%d", &hull.n)!=EOF){
        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.construct();
        pt bc1 = hull.BaryCenter();
        scanf("%d",&hull2.n);
        for (int i = 0; i < hull2.n; i++)
            scanf("%lf%lf%lf", &hull2.P[i].x, &hull2.P[i].y, &hull2.P[i].z);
        hull2.construct();
        pt bc2 = hull2.BaryCenter();
        //printf("BARY1: %.2f %.2f %.2f\n",bc1.x,bc1.y,bc1.z);
        //printf("BARY2: %.2f %.2f %.2f\n",bc2.x,bc2.y,bc2.z);
        double dis1 = Mod, dis2 = Mod;
        for (int i = 0; i < hull.cnt; i++)
            dis1 = min(dis1,fabs(hull.ptoplane(bc1,hull.F[i])));
        for (int i = 0; i < hull2.cnt; i++)
            dis2 = min(dis2,fabs(hull2.ptoplane(bc2,hull2.F[i])));
        printf("%.6f\n",dis1+dis2);
    }
    return 0;
}

时间: 2024-10-13 16:19:57

POJ 2225 / ZOJ 1438 / UVA 1438 Asteroids --三维凸包,求多面体重心的相关文章

UVA 11769 All Souls Night 三维凸包求表面积

题目链接:点击打开链接 求给定的3维坐标的凸包的表面积 三维凸包裸题.. #include <stdio.h> #include <string.h> #include <stdlib.h> #include <math.h> #include <iostream> #include <queue> #include <algorithm> using namespace std; #define PR 1e-8 #def

zoj 1453 Surround the Trees(凸包求周长)

链接:http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=453 Time Limit: 2 Seconds      Memory Limit: 65536 KB There are a lot of trees in an area. A peasant wants to buy a rope to surround all these trees. So at first he must know the minimal r

UVA1438 Asteroids(增量法求三维凸包,加权所有三棱锥质量求多面体重心)

https://www.luogu.com.cn/problem/UVA1438 题解建议参考lrj白书. 代码来自牛逼网友. 我只是存个板子. 1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 #include <cstdlib> 5 #include <vector> 6 #include <algorithm> 7 8 using namespace st

HDU 1240——Asteroids!(三维BFS)POJ 2225——Asteroids

普通的三维广搜,需要注意的是输入:列,行,层 #include<iostream> #include<cstdio> #include<cstring> #include<queue> #include<algorithm> #define M 11 using namespace std; int dir[6][3]={{0,1,0},{0,-1,0},{1,0,0},{-1,0,0},{0,0,1},{0,0,-1}};//6个方向 int

POJ 1504,ZOJ 2001,UVA 713, Adding Reversed Numbers,错误,已找到错误

------------------------------------------------------------ 以此题警告自己: 总结, 1.在数组的使用时,一定别忘了初始化 2.在两种情况复制代码时,一定要小心,注意修改变量名,一不留神就会带来不可估量的后果,一定要仔细挨着一个一个变量的修改,别跳着看着哪个变量就改哪一个变量! (这个题目中,就是复制了一下,代码,ca,我找了一下午的错....还好终于找到了,一个字母的错,) -----------------------------

poj 1979 &amp;&amp; zoj 2165 Red and Black

Red and Black Time Limit: 1000MS   Memory Limit: 30000K Total Submissions: 22409   Accepted: 12100 Description There is a rectangular room, covered with square tiles. Each tile is colored either red or black. A man is standing on a black tile. From a

UVA 152-Tree&#39;s a Crowd(暴力求解三维坐标求最短距离)

Tree's a Crowd Time Limit:3000MS     Memory Limit:0KB     64bit IO Format:%lld & %llu Submit Status Description  Tree's a Crowd  Dr William Larch, noted plant psychologist and inventor of the phrase ``Think like a tree--Think Fig'' has invented a new

ZOJ 2711 Regular Words (三维Catalan数)

题目链接:ZOJ 2711 Regular Words (三维Catalan数) 题意:给出一串3*n长度的字符,其中是含有A,B,C,满足以下条件的字符串有多少种. 1.字符串中A,B,C的个数相同. 2.该字符串的前缀中 A,B,C的个数成非递减. 做法一:一个三维的Catalan数. 做法二:三维的DP. AC代码: import java.math.*; import java.util.*; public class Main{ public static void main(Stri

POJ 3026 Borg Maze &amp; UVA 10307 Killing Aliens in Borg Maze(BFS,最小生成树)

http://poj.org/problem?id=3026 http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=1248 Borg Maze Time Limit: 1000MS   Memory Limit: 65536K Total Submissions: 8498   Accepted: 2862 Description The Bor