博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
三维计算几何之三维凸包 增量法
阅读量:3968 次
发布时间:2019-05-24

本文共 3333 字,大约阅读时间需要 11 分钟。

三维凸包

三维凸包就是将凸包放在三维中求。在三维空间中有一堆点,现求一个多面体将所有点全部包住的最小凸多面体。这里可以类比一下。

这个凸多面体有一个性质:他是包住所有点的多面体中表面积最小的。

三维凸包的求法

三维凸包的求法可以类比二维凸包,但三维凸包的核心还是暴力,暴力的去求。这里主要用一个增量法的思想。

这里要清楚组成三维凸包的不是二维凸包中的一条条线段,在这是一个个平面,由这些平面组成的凸多面体就是三维凸包。
和二维一样,我们先确定一个平面,然后不断地往里加点,每次加一个点,更新一下凸包,直到所有点都加进来了,整个的凸包也就求完了。这就是增量法的思想,不断地一个一个加点。

然后问题就是如何更新凸包呢,使得加入一个点后,更新为一个新的多面体,是i+1个点的凸包。

eg
如上图:每当新加一个点pr时,将这个点想成一个光源(太阳),向四周散发光线(如图二),照到原先的凸多面体CH上,所有能被照到的平面是要删去的面,所有照不到的面是要留下来的面。然后看所有边,每条边都在两个平面上(两平面的交线),当光线照到他所在平面时,代表该边被照到一次,然后我们找出所有被照到一次的边,说明他们是能照到和不能照到的分界线,所以我们新加一个该边和新加入的pr点所构成的一个平面(如图三),就是新凸多面体。

大体的步骤就是如上所说的,但还有很多细节需要注意:存每个平面时都只存三角形平面,即存三个顶点,但若存在n边形的平面怎么办,我们可以把他划分成n-2个三角形,然后存下来,因为任何的多边形都可以通过三角剖分得到多个三角形。

其次如何判断一个点能照到的平面呢?我们可以判断这个点在平面上方还是下方(通过和法向量的点积判断,具体在),在上面说明能照到,在下面说明照不到,所以这里的平面是有方向的,因此我们存平面时一定要严格按照一个方向存三个点,这里统一规定逆时针存。
然后就是判断每条边被照到了几次,我们加个bool二维数组,g[a][b]为1的话说明a到b的边被照到一次,所以当g[a][b]为1,g[b][a]为0时说明ab这条边只被照到了一次。(因为都是逆时针存点,所以这个面是ab边时,邻面就是ba边)
还有就是如何只存三角形平面呢?我们给每个点都加一个微小扰动,就是让他的坐标发生一个很细微很细微的变化,使得不存在四点共面的情况,然后任意三个点就可以构成一个唯一平面了。(注意加了微小扰动之后就不能用a和b的绝对值差eps内就算相等了,必须严格大于等于小于判断)

下面具体可以看代码模板了解详细步骤


代码模板

#include
#include
#include
using namespace std;const int N = 210;const double eps = 1e-12; //精度要求高一点int n,m; //n是总点数,m是凸包中平面的数量bool g[N][N]; //g用来判断一条边被照到几次double rand_eps() //用rand函数来生成一个非常小的随机数{
return ((double)rand() / RAND_MAX - 0.5) * eps; //用rand生成一个-0.5到0.5之间的数,再乘eps,就得到了一个非常小的随机数}struct Point{
//定义点的结构体 double x,y,z; //xyz三个坐标 void shake() //微小扰动,给每个坐标都加一个极小的随机数 {
x += rand_eps(), y += rand_eps(), z += rand_eps(); } Point operator-(Point t) //重载一下减号运算符 {
return {
x-t.x, y-t.y, z-t.z}; } Point operator*(Point t) //叉乘 {
return {
y*t.z - t.y*z, t.x*z - x*t.z, x*t.y - y*t.x}; } double operator&(Point t) //点积 {
return x*t.x + y*t.y + z*t.z; } double len() //求向量的模长,三个坐标也能存向量 {
return sqrt(x*x + y*y + z*z); }}p[N]; //p来存所有点struct Plane{
//定义平面的结构体 int v[3]; //三个顶点 Point norm() //求法向量 {
return (p[v[1]] - p[v[0]]) * (p[v[2]] - p[v[0]]); //平面中两向量的叉积 } bool above(Point t ) //判断一个点是否在平面上方 {
return ((t-p[v[0]]) & norm()) >= 0; //用向量和法向量的点积判断 } double area() //求一平面的面积 {
return norm().len() / 2; //法向量的模长除以2 }}plane[N],tp[N]; //plane存凸包上的平面,tp用来更新凸包,每次凸包上要留的平面和新加的平面都存进tp,要删的平面不存,最后将tp在复制给plane,就实现了凸包的更新void convex(){
plane[m++] = {
0,1,2}; //初始化凸包,随便三个点存入,确定最开始的一个平面,这里取得是前三个点 plane[m++] = {
2,1,0}; //因为不知道第一个平面怎么样是逆时针,所以都存一遍,顺时针存的一会会被删掉 for(int i = 3;i < n;i++) //从第四个点开始循环每个点 {
int cnt = 0; for(int j = 0;j < m;j++) //循环每个平面 {
bool fg = plane[j].above(p[i]); //判断这个点是否在该平面上方 if(!fg) //如果是下方的话,说明照不到 tp[cnt++] = plane[j]; //存进tp数组 for(int k = 0;k < 3;k++) //循环该平面的三条边 g[plane[j].v[k]][plane[j].v[(k+1)%3]] = fg; //ab边照不照得到情况赋值给g[a][b] } for(int j = 0;j < m;j++) //然后就循环每个平面的每条边 {
for(int k = 0;k < 3;k++) {
int a = plane[j].v[k],b = plane[j].v[(k+1)%3]; if(g[a][b] && !g[b][a]) //判断该边是否被照到了一次,即是否是交界线的边 tp[cnt++] = {
a,b,i}; //若是,加新平面abi,ab一定是逆时针的,i在后面 } } m = cnt; //将tp再赋值给plane for(int j = 0;j < m;j++) plane[j] = tp[j]; }}int main(){
cin >> n; for(int i = 0;i < n;i++) {
cin >> p[i].x >> p[i].y >> p[i].z; //输入n个点 p[i].shake(); //微小扰动 } convex(); //求三维凸包 double ans = 0; //求面积和 for(int i = 0;i < m;i++) //循环最终凸包上的m个平面 ans += plane[i].area(); //将平面的面积加和 printf("%.6lf\n",ans); //输出答案 return 0;}

转载地址:http://abbki.baihongyu.com/

你可能感兴趣的文章
Java 集合框架
查看>>
Weblogic 精萃
查看>>
Servlet 精萃
查看>>
XStream 精萃
查看>>
XStream 环境设置
查看>>
Git 分支
查看>>
Git 冲突
查看>>
Git Merging vs. Rebasing
查看>>
libreoffice/openoffice c/c++转换office格式为pdf
查看>>
Tomcat 7.0 64位免安装解压版 安装及配置
查看>>
Android 网络编程 初级入门(一)
查看>>
No enclosing instance of type Demo06 is accessible.
查看>>
计算机发展中的两大“杀手”
查看>>
MDK5(Keil for ARM) 工程建立时遇到的问题集锦
查看>>
Ubuntu下安装GTK+及Glade开发C应用界面
查看>>
assertion 'GTK_IS_WIDGET (widget)' failed的解决办法
查看>>
Ubuntu登录管理员账户时,输入密码后一直在登录界面循环
查看>>
Linux下的定时器以及POSIX定时器:timer_settime()
查看>>
POSIX定时器timer_create()以及线程中的gettid() 和pthread_self()
查看>>
c /c++中日期和时间的获取:strftime()函数
查看>>