三維凸包即解決空間中n個點構成的最小空間幾何體的相關問題
選擇題目：洛谷P4724
我所學習的方法是增量法 （與最小覆寫圓類似） 首先選出四個不共面的點為基礎凸包
不斷加點 判斷這個點是否在這個凸包內 不在就形成新的凸包
具體的看代碼理解
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MAXN=2050;
const double eps=1e-8;
struct Point3{
//三維的點
double x,y,z;
Point3(){}
Point3(double x,double y,double z):x(x),y(y),z(z){}///多載三維 點的計算
Point3 operator + (Point3 B){return Point3(x+B.x,y+B.y,z+B.z);}
Point3 operator - (Point3 B){return Point3(x-B.x,y-B.y,z-B.z);}
Point3 operator * (double k){return Point3(x*k,y*k,z*k);}
Point3 operator / (double k){return Point3(x/k,y/k,z/k);}
};
typedef Point3 Vector3;
double Dot(Vector3 A,Vector3 B){return A.x*B.x+A.y*B.y+A.z*B.z;}///三維點積
Point3 Cross(Vector3 A,Vector3 B){return Point3(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x);}///三維叉積
double Len(Vector3 A){return sqrt(Dot(A,A));}
///向量的長度
double Area2(Point3 A,Point3 B,Point3 C){return Len(Cross(B-A,C-A));}
///四面體有向體積*6
double volume4(Point3 A,Point3 B,Point3 C,Point3 D){return Dot(Cross(B-A,C-A),D-A);}
struct CH3D{
struct face{
ll a,b,c;
///凸包的一個面上的3個點的編號
bool ok;
///該面是否在最終的凸包上
};
ll n;
///初始頂點數
Point3 P[MAXN];
///初始頂點
ll num;
///凸包表面的三角形個數
face F[8*MAXN];
///凸包表面的三角形
ll g[MAXN][MAXN];/// g[x][y]表示點x到點y屬于哪個面
///點在面的同向
double dblcmp(Point3 &p,face &f){
Point3 m=P[f.b]-P[f.a];
Point3 n=P[f.c]-P[f.a];
Point3 t=p-P[f.a];
return Dot(Cross(m,n),t);
}
void deal(ll p,ll a,ll b){
ll f=g[a][b];
///搜索與該邊相鄰的另一個平面
face add;
if(F[f].ok){
if(dblcmp(P[p],F[f])>eps){
///如果如果從p點能看到該面f,則繼續深度搜索f的三條邊
///以更新新的凸面（維護）
F[f].ok=false;
deal(p,F[f].b,F[f].a);
deal(p,F[f].c,F[f].b);
deal(p,F[f].a,F[f].c);
}
else{
///如果如果從p點不能看到f面，則p點和a、b點組成一個三角形
add.a=b;
add.b=a;
add.c=p;
add.ok=true;
g[p][b]=g[a][p]=g[b][a]=num;
F[num++]=add;
}
}
}
void dfs(ll p,ll now){
///維護凸包，如果點p在凸包外則更新凸包
F[now].ok=false;
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 same(ll s,ll t){
///判斷兩個面是否為同一面
Point3 &a=P[F[s].a];
Point3 &b=P[F[s].b];
Point3 &c=P[F[s].c];
return fabs(volume4(a,b,c,P[F[t].a]))<eps &&
fabs(volume4(a,b,c,P[F[t].b]))<eps &&
fabs(volume4(a,b,c,P[F[t].c]))<eps ;
}
void create(){
///構建三維凸包
ll i,j,tmp;
face add;
num=0;
if(n<4)return;
bool flag=true;
for(i=1;i<n;i++){
///使前兩個點不共點
if(Len(P[0]-P[i])>eps){
swap(P[1],P[i]);
flag=false;
break;
}
}
if(flag)return;
flag=true;
for(i=2;i<n;i++){
///使前三個點不共線
if(Len(Cross(P[0]-P[1],P[1]-P[i]))>eps){
swap(P[2],P[i]);
flag=false;
break;
}
}
if(flag)return;
flag=true;
for(ll i=3;i<n;i++){
///使前四個點不共面
if(fabs(Dot(Cross(P[0]-P[1],P[1]-P[2]),P[0]-P[i]))>eps){
swap(P[3],P[i]);
flag=false;
break;
}
}
if(flag)return;
for(i=0;i<4;i++){
///構建初始四面體
///四個點為：
///p[0],p[1],p[2],p[3]
add.a=(i+1)%4;
add.b=(i+2)%4;
add.c=(i+3)%4;
add.ok=true;
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]=num;
///逆向的有向邊保存
F[num++]=add;
}
for(i=4;i<n;i++){
///構建&更新凸包
for(j=0;j<num;j++){
///判斷點是否在當前三維凸包內，i表示當前點，j表示當前面
if(F[j].ok && dblcmp(P[i],F[j])>eps){
///對當前凸包面進行判斷，看點能否看到這個面
dfs(i,j);
///點能看到當前面，更新凸包的面
break;
}
}
}
tmp=num;
for(i=num=0;i<tmp;i++){
if(F[i].ok)F[num++]=F[i];
}
}
///凸包的表面積
double area(){
double res=0;
for(ll i=0;i<num;i++)res+=Area2(P[F[i].a],P[F[i].b],P[F[i].c]);
return res/2.0;
}
///體積
double volume(){
double res=0;
Point3 tmp(0,0,0);
for(ll i=0;i<num;i++)res+=volume4(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
return fabs(res/6.0);
}
///表面三角形個數
ll triangle(){
return num;
}
///表面多邊形個數
ll polygon(){
ll i,j,res,flag;
for(i=res=0;i<num;i++){
flag=true;
for(j=0;j<i;j++){
if(same(i,j)){
flag=false;
break;
}
}
res+=flag;
}
return res;
}
};
CH3D hull;
int main(){
scanf("%lld",&hull.n);
for(ll i=0;i<hull.n;i++)scanf("%lf %lf %lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
hull.create();
printf("%.3lf\n",hull.area());
///問啥輸啥就完了
return 0;
}
看到最后 來個彩蛋 推薦寫博客正在聽的輕音樂 ：Silent Night
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/55356.html
標籤：其他
上一篇：Codeforces Round #637 (Div. 2) A~D
下一篇：資料結構（C語言版）---查找