计算 3D 网格边界框对角线长度的简单方法是什么? [英] What's the simple way to compute the diagonal length of a 3D mesh bounding box?
问题描述
我想计算 3D 网格边界框的对角线长度.使用 C++，我迭代顶点并搜索 X 坐标的 (min,max)、Y 坐标的 (min,max) 和 Z 坐标的 (min,max).但是我不知道如何利用这些获得的最小值/最大值来计算边界框的对角线长度.有什么帮助吗?
为简单起见，让我们考虑将 n 个 3D 点(点云)作为输入(而不是网格)的列表，这对于多边形网格.
网格的对角线"就是网格中两个最远点之间的线.这很容易用简单的 O(n^2) 蛮力搜索(2 个嵌套的 for 循环记住最远的点)来计算.还有一些更快的方法可以利用点的排序.这里是蛮力示例:
line pointcloud::diagonal(){int i,j;线 l,ll;l=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0));//空行for (i=0;i<pnt.num-1;i++)//O(n^2) 搜索所有点对对于 (j=i+1;j有关 line 和 pointcloud 类实现的更多信息，请阅读下面的链接和 OBB 的源代码.
但是，从评论中我感觉到您需要 3D OBB(定向边界框)而不仅仅是对角线.您现在拥有的只是 AABB(轴对齐的边界框)，它不会为您提供网格对角线(除非它的幸运方向与 AABB 对角线匹配).>
注意 AABB 和 OBB 对角线与网格对角线不一样！！！
有很多方法可以计算OBB，从蛮力(~O(n^6))到使用特征向量、凸包等更快...
我设法将我的
求最小边界体积
只需生成某个坐标系(覆盖半球)的所有 m 旋转.你不需要覆盖整个球体，因为另一半只是否定......
现在通过获取沿两个方向的 3 个轴的距离并计算形成的盒子的体积并记住最小的一个(轴、距离和体积)来计算每个体积.由于舍入和非线性问题，在 cube_map 中可能有单元化的数据导致 volume = 0(如果 cube_map 在开始时被清除为零)所以忽略这样的体积.
在此之后，您应该获得OBB 近似值.以下是几个旋转位置的 OBB 预览:
它有点跳跃，因为对于这种对称形状，有效的OBB数量是无限的，并且在不同的旋转中可以在搜索中首先找到不同的.
提高准确性
只需搜索附近的几个旋转找到OBB 近似值并记住最小的一个.这可以递归地完成.但是我懒得去实现这个，因为 OBB 结果的当前状态对我来说已经足够了.
这里是 C++/GL 源代码(其余的可以在上面的链接中找到):
//---------------------------------------------------------------------类点云{民众://配置文件列表点;点云(){}pointcloud(pointcloud& a) { *this=a;}〜点云(){}pointcloud* operator = (const pointcloud *a) { *this=*a;返回这个；}//pointcloud* operator = (const pointcloud &a) { ...copy... return this;}void reset(){ pnt.num=0;}void add(vec3 p){ pnt.add(p);}void add(point p){ pnt.add(p.p0);}无效计算(){}；无效绘制(){glBegin(GL_POINTS);for (int i=0;i类cube_map{民众:int n,nn,sz;浮动 fn2;T图[6*N*N]；立方体地图(){ n = N;nn=N*N；sz=6*nn;fn2=0.5*float(n);}立方体地图(立方体地图& a){ *this=a;}~cube_map() {}cube_map* 运算符 = (const cube_map *a) { *this=*a;返回这个；}//cube_map* operator = (const cube_map &a) { ...copy... return this;}vec3 ix2dir(int ix){浮动 x,y;vec3 dir=vec3(0.0,0.0,0.0);if ((ix<0)||(ix>=sz)) 返回目录；x=ix%n;ix/=n;x/=fn2;X - ;y=ix%n;ix/=n;y/=fn2;y--;如果 (ix==0){ dir.y=x;dir.z=y;目录.x=-1.0;}如果 (ix==1){ dir.y=x;dir.z=y;目录.x=+1.0;}如果 (ix==2){ dir.x=x;dir.z=y;目录.y=-1.0;}如果 (ix==3){ dir.x=x;dir.z=y;目录.y=+1.0;}如果 (ix==4){ dir.x=x;目录.y=y;dir.z=-1.0;}如果 (ix==5){ dir.x=x;目录.y=y;dir.z=+1.0;}返回归一化(目录)；}int dir2ix(vec3 目录){int ix=0,x=0,y=0;浮动 a=0.0,b;b=fabs(dir.x);如果 (a<b){ a=b;如果 (dir.x<0) ix=0;否则 ix=1;}b=fabs(dir.y);如果 (a<b){ a=b;如果 (dir.y<0) ix=2;否则 ix=3;}b=fabs(dir.z);如果 (a<b){ a=b;如果 (dir.z<0) ix=4;否则 ix=5;}目录/=a;dir+=vec3(1.0,1.0,1.0);目录*=fn2;if (ix==0){ x=dir.y;y=dir.z;}if (ix==1){ x=dir.y;y=dir.z;}if (ix==2){ x=dir.x;y=dir.z;}if (ix==3){ x=dir.x;y=dir.z;}if (ix==4){ x=dir.x;y=dir.y;}if (ix==5){ x=dir.x;y=dir.y;}ix=(ix*nn)+(y*n)+(x)；如果 ((ix<0)||(ix>=sz)) ix=0;返回 ix；}void set(vec3 dir,T &a){ map[dir2ix(dir)]=a;}T get(vec3 dir){ return map[dir2ix(dir)];}void clear(T &a){ for (int i=0;i<sz;i++) map[i]=a;}};//---------------------------------------------------------------------------class OBB//定向边界框{民众://计算vec3 p0;//中央vec3 u,v,w;//基半向量(p0 原点)OBB() {}OBB(OBB& a) { *this=a;}~OBB() {}OBB* 运算符 = (const OBB *a) { *this=*a;返回这个；}//OBB* operator = (const OBB &a) { ...copy... return this;}无效计算(点云和pcl){const int N=24;int i,j,k,na=6*N,nb=2*N;cube_map地图;垫4米，马;vec3 o,p,q,pp0;国际a,b;浮动 da,db,d,dd,e,ee,V,VV;p0=vec3(0.0,0.0,0.0);u=vec3(0.0,0.0,0.0);v=vec3(0.0,0.0,0.0);w=vec3(0.0,0.0,0.0);如果 (pcl.pnt.num<=0) 返回；//初始化常量和东西da=2.0*M_PI/float(na);db= M_PI/float(nb-1);//计算平均点对于 (j=0;j未旋转坐标系(1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0);对于 ( a=0;a<na;a+=2,ma=lrotz(ma,da))对于 (m=lroty(ma,float(-0.5*M_PI)),b=0;b<nb;b++,m=lroty(m,db)){//获取每个 m 方向的 OBBp.x=map.get(-m[0].xyz);q.x=map.get(+m[0].xyz);p.y=map.get(-m[1].xyz);q.y=map.get(+m[1].xyz);p.z=map.get(-m[2].xyz);q.z=map.get(+m[2].xyz);o=p+q；VV=fabs(o.x*o.y*o.z);如果 ((V>VV)&&(VV>1e-6)){V=VV；u=m[0].xyz;v=m[1].xyz;w=m[2].xyz;o*=0.5；pp0=p0+(u*(o.x-p.x))+(v*(o.y-p.y))+(w*(o.z-p.z));u*=o.x;v*=o.y;w*=o.z;}}p0=pp0;}无效绘制(){const vec3 p[8]={p0-u-v-w,p0+u-v-w,p0+u+v-w,p0-u+v-w,p0-u-v+w,p0+u-v+w,p0+u+v+w,p0-u+v+w,};const int ix[24]={0,1,1,2,2,3,3,0,4,5,5,6,6,7,7,4,0,4,1,5,2,6,3,7,};glBegin(GL_LINES);for (int i=0;i<24;i++) glVertex3fv(p[ix[i]].dat);glEnd();}};//--------------------------------------------------------------------------- 希望我没有忘记复制一些东西......我想尽可能地保持代码简单，所以它不是很优化，并且有很大的改进空间.使用的数学是基于 GLSL 的，因此您可以使用 GLM.我使用了我自己的库，如果需要，可以在上面的链接中找到 vec(但需要生成其 ~220KByte 的代码)但完全匹配 GLSL 和 GLM，所以你可以使用它.然而，mat4 使用了一些 GLM 中不存在的这种格式的函数，以防万一:
template _mat4 类{民众:_vec4颜色 [4];//列！！！_mat4(T a00,T a01,T a02,T a03,T a04,T a05,T a06,T a07,T a08,T a09,T a10,T a11,T a12,T a13,T a14,T a15){col[0]=vec4(a00,a01,a02,a03);//x 轴col[1]=vec4(a04,a05,a06,a07);//y 轴col[2]=vec4(a08,a09,a10,a11);//z轴col[3]=vec4(a12,a13,a14,a15);//起源}_mat4(){col[0]=vec4(1,0,0,0);col[1]=vec4(0,1,0,0);col[2]=vec4(0,0,1,0);col[3]=vec4(0,0,0,1);}_mat4(const _mat4& a) { *this=a;}~_mat4() {}//运算符(矩阵数学)_mat4* operator = (const _mat4 &a) { for (int i=0;i<4;i++) col[i]=a.col[i];返回这个；}//=a[][]_vec4<T&运算符 [](const int i){ 返回 col[i];}//a[i]_mat4运算符 * (_mat4<T>&m)//=a[][]*m[][]{_mat4q;输入 i,j,k;对于 (i=0;i<4;i++)对于 (j=0;j<4;j++)对于 (q.col[i][j]=0,k=0;k<4;k++)q.col[i][j]+=col[k][j]*m.col[i][k];返回 q；}_mat4运算符 * (_vec4<T>&v)//=a[][]*v[]{_vec4q;int i,j;对于 (i=0;i<4;i++)对于 (q.dat[i]=0,j=0;j<4;j++)q.dat[i]+=col[i][j]*v.dar[j];返回 q；}_mat4运算符 * (T &c)//=a[][]*c{_mat4q;int i,j;对于 (i=0;i<4;i++)对于 (j=0;j<4;j++)q.dat[i]=col[i][j]*c;返回 q；}_mat4运算符/(T &c)//=a[][]/c{_mat4q;int i,j;对于 (i=0;i<4;i++)对于 (j=0;j<4;j++)q.dat[i]=divide(col[i][j],c);返回 q；}_mat4运算符 *=(_mat4<T>&m){ this[0]=this[0]*m;返回 *this;};_mat4运算符 *=(_vec4<T>&v){ this[0]=this[0]*v;返回 *this;};_mat4运算符 *=(const T &c){ this[0]=this[0]*c;返回 *this;};_mat4运算符/=(const T &c){ this[0]=this[0]/c;返回 *this;};//成员void get(T *a){输入 i,j,k;对于 (k=0,i=0;i<4;i++)对于 (j=0;j<4;j++,k++)a[k]=col[i].dat[j];}空集(T * a){输入 i,j,k;对于 (k=0,i=0;i<4;i++)对于 (j=0;j<4;j++,k++)col[i].dat[j]=a[k];}};//---------------------------------------------------------------------------模板_mat4转置(const _mat4<T> &m){_mat4q;int i,j;对于 (i=0;i<4;i++)对于 (j=0;j<4;j++)q.col[i][j]=m.col[j][i];返回 q；}//---------------------------------------------------------------------------模板_mat4逆(_mat4<T>&m){Tp[3];_mat4q;T x,y,z;int i,j;//转置旋转for (i=0;i<3;i++) for (j=0;j<3;j++) q.col[i][j]=m.col[j][i];//复制投影对于 (i=0;i<4;i++) q.col[i][3]=m.col[i][3];//转换原点:new_pos = - new_rotation_matrix * old_posfor (i=0;i<3;i++) for (p[i]=0,j=0;j<3;j++) p[i]+=q.col[j][i]*m.col[3][j];对于 (i=0;i<3;i++) q.col[3][i]=-p[i];返回 q；}//---------------------------------------------------------------------------模板_mat4lrotx(_mat4<T>&m,Tang){Tc=cos(ang),s=sin(ang);_mat4r=mat4(1, 0, 0, 0,0, c, s, 0,0,-s, c, 0,0, 0, 0, 1);r=m*r；返回 r;};//---------------------------------------------------------------------------模板_mat4lroty(_mat4<T>&m,Tang){Tc=cos(ang),s=sin(ang);_mat4r=mat4(c, 0,-s, 0,0, 1, 0, 0,s, 0, c, 0,0, 0, 0, 1);r=m*r；返回 r;};//---------------------------------------------------------------------------模板_mat4lrotz(_mat4<T>&m,Tang){Tc=cos(ang),s=sin(ang);_mat4r=mat4(c, s, 0, 0,-s, c, 0, 0,0, 0, 1, 0,0, 0, 0, 1);r=m*r；返回 r;};//---------------------------------------------------------------------------模板_mat4旋转(_mat4 &m,T ang,_vec3 p0,_vec3 dp){国际我;Tc=cos(ang),s=sin(ang);_vec3x,y,z;_mat4a,_a,r=mat4(1, 0, 0, 0,0, c, s, 0,0,-s, c, 0,0, 0, 0, 1);//基向量x=归一化(dp)；//旋转轴y=_vec3 T (1,0,0)；//任何不平行于 x 的向量如果(fabs(dot(x,y))＞0.75)y＝_vec3 T(0,1,0)；z=cross(x,y);//z 垂直于 x,yy=cross(z,x);//y 垂直于 x,zy=归一化(y)；z=归一化(z)；//输入矩阵对于 (i=0;i<3;i++){a[0][i]=x[i];a[1][i]=y[i];a[2][i]=z[i];a[3][i]=p0[i];a[i][3]=0;} a[3][3]=1;_a=逆(a)；r=m*a*r*_a;返回 r;};//---------------------------------------------------------------------------模板_mat4grotx(_mat4<T>&m,Tang){ return inverse(lrotx(inverse(m),ang));};模板_mat4groty(_mat4<T> &m,T ang){ return inverse(lroty(inverse(m),ang));};模板_mat4grotz(_mat4<T>&m,Tang){ return inverse(lrotz(inverse(m),ang));};//---------------------------------------------------------------------------typedef _mat4垫4;typedef _mat4dmat4;typedef _mat4bmat4;typedef _mat4imat4;typedef _mat4umat4;//---------------------------------------------------------------------------mat4 GLSL_math_test4x4;//--------------------------------------------------------------------------- 为了理解它或写你自己的，我建议看:
最后我还使用了我的动态列表模板:
List 等同于 double xxx[];xxx.add(5); 将 5 添加到列表末尾xxx[7] 访问数组元素(安全)xxx.dat[7] 访问数组元素(不安全但快速的直接访问)xxx.num 是数组实际使用的大小xxx.reset() 清空数组并设置xxx.num=0xxx.allocate(100) 为 100 个项目预分配空间
现在OBB的结果
它只是由其中心 p0 和半向量 u,v,w 描述的框.所以要获得点云的OBB PCL 只需计算:
OBB obb;点云PCL；PCL.reset();PCL.add(...);//这里将点输入 PCLobb.compute(PCL);仅此而已.
I would like to compute the diagonal length of the bounding box of a 3D mesh. Using C++, I iterate of the vertices and search for the (min,max) of X coordinates, (min,max) of Y coordinates and (min,max) of Z coordinates. However I don't know how can I exploit these obtained min/max for the computation of the diagonal length of the bounding box. Any help please?
For simplicity let us consider a list of n 3D points (point cloud) as input (instead of mesh) which is enough for polygonal meshes.
The "diagonal" of the mesh is just line between 2 most distant points in the mesh. That is easily computable with trivial O(n^2) brute force search (2 nested for loops remembering most distant points). There are also faster methods that exploit ordering of points. Here the brute force example:
line pointcloud::diagonal()
{
int i,j;
line l,ll;
l=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0)); // empty line
for (i=0;i<pnt.num-1;i++) // O(n^2) search through all point pairs
for (j=i+1;j<pnt.num-1;j++)
{
ll=line(pnt.dat[i],pnt.dat[j]); // prepare line
if (l.l<ll.l) l=ll; // compare sizes and remember the longer one
}
return l;
}
For more info about line and pointcloud classes implementation read the links below and source code for the OBB.
However from the comments I got the feeling you need 3D OBB (oriented bounding box) instead of just diagonal. What you have right now is just AABB (axis aligned bounding box) which will not give you the mesh diagonal (unless its in lucky orientation that matches AABB diagonal).
Beware both AABB and OBB diagonal is not the same as mesh diagonal !!!
There are many methods to compute OBB from brute force (~O(n^6)) to faster using eigen vectors, convex hull etc...
I managed to port my 2D OBB approximation into 3D.
The idea is the same. Store max distances in "all" (m) possible directions/angles (covering full sphere instead of circle in 2D) reducing data from n to m. And then just search the computed data for minimal bounding volume (instead of area in 2D).
I used my Cone to box collision for testing and as a start point.
The algo:
compute pivot point
p0it must be inside point of the OBB. usually center of AABB or avg point is enough for this.
compute distances in each possible direction
there is infinite number of possible directions so we need to limit this to
m. the bigger themthe slower computation but more accurate. In order to store and obtain these values fast I usedcube_map.Its a 2D texture covering the surface of unit cube (6 x square sides) and its addressed by direction vector instead of texture coordinates.
I implemented 2 functions that convert between
indexin texture data (stored as 1D array) anddirectionvector. For more info seecube_mapin the example...The distance
dof pointpfromp0in some directiondiris computed like this:d = dot( p-p0 , dir )so generate
mpossible directions, and for each compute distance for all points in your source list of point and remember the biggest one which is then stored tocube_mapfor latter. This isO(m*n)Here example of stored distances for one frame (content of cube_map):
find minimal bounding volume
Simply generate all
mrotations of some coordinate system (covering half sphere). You do not need to cover full sphere because the other half is just negation...Now for each compute volume by getting the distances along its 3 axises in both directions and computing the volume of formed box and remember the smallest one (axises, distances and volume). There is possibility of having unitialized data in the
cube_mapwhich results involume = 0(if cube_map was cleared to zero at start) due rounding and nonlinearity problems so ignore such just volumes.After this you should have your OBB aproximation. Here preview of OBB for few rotated positions:
Its a bit jumpy because for such symmetric shape there are infinite number of valid OBBs and in different rotations different one can be found first in search.
improve accuracy
Simply search few rotations nearby found OBB aproximation and remember the smallest one. This can be done recursively. However I am too lazy to implement this as current state of OBB result is enough for me.
Here C++/GL source (the rest can be found in the link above):
//---------------------------------------------------------------------------
class pointcloud
{
public:
// cfg
List<vec3> pnt;
pointcloud() {}
pointcloud(pointcloud& a) { *this=a; }
~pointcloud() {}
pointcloud* operator = (const pointcloud *a) { *this=*a; return this; }
//pointcloud* operator = (const pointcloud &a) { ...copy... return this; }
void reset(){ pnt.num=0; }
void add(vec3 p){ pnt.add(p); }
void add(point p){ pnt.add(p.p0); }
void compute(){};
void draw()
{
glBegin(GL_POINTS);
for (int i=0;i<pnt.num;i++) glVertex3fv(pnt.dat[i].dat);
glEnd();
}
};
//---------------------------------------------------------------------------
template<class T,int N> class cube_map
{
public:
int n,nn,sz;
float fn2;
T map[6*N*N];
cube_map() { n=N; nn=N*N; sz=6*nn; fn2=0.5*float(n); }
cube_map(cube_map& a) { *this=a; }
~cube_map() {}
cube_map* operator = (const cube_map *a) { *this=*a; return this; }
//cube_map* operator = (const cube_map &a) { ...copy... return this; }
vec3 ix2dir(int ix)
{
float x,y;
vec3 dir=vec3(0.0,0.0,0.0);
if ((ix<0)||(ix>=sz)) return dir;
x=ix%n; ix/=n; x/=fn2; x--;
y=ix%n; ix/=n; y/=fn2; y--;
if (ix==0){ dir.y=x; dir.z=y; dir.x=-1.0; }
if (ix==1){ dir.y=x; dir.z=y; dir.x=+1.0; }
if (ix==2){ dir.x=x; dir.z=y; dir.y=-1.0; }
if (ix==3){ dir.x=x; dir.z=y; dir.y=+1.0; }
if (ix==4){ dir.x=x; dir.y=y; dir.z=-1.0; }
if (ix==5){ dir.x=x; dir.y=y; dir.z=+1.0; }
return normalize(dir);
}
int dir2ix(vec3 dir)
{
int ix=0,x=0,y=0;
float a=0.0,b;
b=fabs(dir.x); if (a<b){ a=b; if (dir.x<0) ix=0; else ix=1; }
b=fabs(dir.y); if (a<b){ a=b; if (dir.y<0) ix=2; else ix=3; }
b=fabs(dir.z); if (a<b){ a=b; if (dir.z<0) ix=4; else ix=5; }
dir/=a;
dir+=vec3(1.0,1.0,1.0);
dir*=fn2;
if (ix==0){ x=dir.y; y=dir.z; }
if (ix==1){ x=dir.y; y=dir.z; }
if (ix==2){ x=dir.x; y=dir.z; }
if (ix==3){ x=dir.x; y=dir.z; }
if (ix==4){ x=dir.x; y=dir.y; }
if (ix==5){ x=dir.x; y=dir.y; }
ix=(ix*nn)+(y*n)+(x);
if ((ix<0)||(ix>=sz)) ix=0;
return ix;
}
void set(vec3 dir,T &a){ map[dir2ix(dir)]=a; }
T get(vec3 dir ){ return map[dir2ix(dir)]; }
void clear(T &a){ for (int i=0;i<sz;i++) map[i]=a; }
};
//---------------------------------------------------------------------------
class OBB // Oriented Bounding Box
{
public:
// computed
vec3 p0; // center
vec3 u,v,w; // basis half vectors (p0 origin)
OBB() {}
OBB(OBB& a) { *this=a; }
~OBB() {}
OBB* operator = (const OBB *a) { *this=*a; return this; }
//OBB* operator = (const OBB &a) { ...copy... return this; }
void compute(pointcloud &pcl)
{
const int N=24;
int i,j,k,na=6*N,nb=2*N;
cube_map<float,N> map;
mat4 m,ma;
vec3 o,p,q,pp0;
int a,b;
float da,db,d,dd,e,ee,V,VV;
p0=vec3(0.0,0.0,0.0);
u=vec3(0.0,0.0,0.0);
v=vec3(0.0,0.0,0.0);
w=vec3(0.0,0.0,0.0);
if (pcl.pnt.num<=0) return;
// init constants and stuff
da=2.0*M_PI/float(na );
db= M_PI/float(nb-1);
// compute avg point
for (j=0;j<pcl.pnt.num;j++) p0+=pcl.pnt.dat[j];
p0/=pcl.pnt.num;
// [compute perpendicular distances]
// fill whole surface of cubemap
for (map.clear(0.0),i=0;i<map.sz;i++)
{
// cube map index to 3D direction
p=map.ix2dir(i);
// compute max distance from p0 in direction p
d=dot(pcl.pnt.dat[0]-p0,p);
for (j=1;j<pcl.pnt.num;j++)
{
dd=dot(pcl.pnt.dat[j]-p0,p);
if (d<dd) d=dd;
}
// store it in cube map for latter
map.map[i]=d;
}
// [pick the smallest volume OBB combination]
V=1e300; pp0=p0;
// try half of "all" rotations (the other one is just negation)
ma=mat4 // unit matrix -> unrotated coordinate system
(
1.0,0.0,0.0,0.0,
0.0,1.0,0.0,0.0,
0.0,0.0,1.0,0.0,
0.0,0.0,0.0,1.0
);
for ( a=0;a<na;a+=2,ma=lrotz(ma,da))
for (m=lroty(ma,float(-0.5*M_PI)),b=0;b<nb;b++,m=lroty(m,db))
{
// get OBB per orientation of m
p.x=map.get(-m[0].xyz);
q.x=map.get(+m[0].xyz);
p.y=map.get(-m[1].xyz);
q.y=map.get(+m[1].xyz);
p.z=map.get(-m[2].xyz);
q.z=map.get(+m[2].xyz);
o=p+q;
VV=fabs(o.x*o.y*o.z);
if ((V>VV)&&(VV>1e-6))
{
V=VV;
u=m[0].xyz;
v=m[1].xyz;
w=m[2].xyz;
o*=0.5;
pp0=p0+(u*(o.x-p.x))+(v*(o.y-p.y))+(w*(o.z-p.z));
u*=o.x;
v*=o.y;
w*=o.z;
}
}
p0=pp0;
}
void draw()
{
const vec3 p[8]=
{
p0-u-v-w,
p0+u-v-w,
p0+u+v-w,
p0-u+v-w,
p0-u-v+w,
p0+u-v+w,
p0+u+v+w,
p0-u+v+w,
};
const int ix[24]=
{
0,1,1,2,2,3,3,0,
4,5,5,6,6,7,7,4,
0,4,1,5,2,6,3,7,
};
glBegin(GL_LINES);
for (int i=0;i<24;i++) glVertex3fv(p[ix[i]].dat);
glEnd();
}
};
//---------------------------------------------------------------------------
Hope I did not forget to copy something... I wanted to keep the code as simple as I could so its not very optimized and there is a lot of room for improvement. The math used is GLSL based so you can use GLM. I used my own libs for that the vec can be found in the links above if needed (but need to be generated as its ~220KByte of code) but is matches GLSL and GLM exactly so you cn use that. The mat4 however use some functions that are not present in GLM in such format so just in case:
template <class T> class _mat4
{
public:
_vec4<T> col[4]; // columns!!!
_mat4(T a00,T a01,T a02,T a03,T a04,T a05,T a06,T a07,T a08,T a09,T a10,T a11,T a12,T a13,T a14,T a15)
{
col[0]=vec4(a00,a01,a02,a03); // x axis
col[1]=vec4(a04,a05,a06,a07); // y axis
col[2]=vec4(a08,a09,a10,a11); // z axis
col[3]=vec4(a12,a13,a14,a15); // origin
}
_mat4()
{
col[0]=vec4(1,0,0,0);
col[1]=vec4(0,1,0,0);
col[2]=vec4(0,0,1,0);
col[3]=vec4(0,0,0,1);
}
_mat4(const _mat4& a) { *this=a; }
~_mat4() {}
// operators (matrix math)
_mat4* operator = (const _mat4 &a) { for (int i=0;i<4;i++) col[i]=a.col[i]; return this; } // =a[][]
_vec4<T>& operator [](const int i){ return col[i]; } // a[i]
_mat4<T> operator * (_mat4<T>&m) // =a[][]*m[][]
{
_mat4<T> q;
int i,j,k;
for (i=0;i<4;i++)
for (j=0;j<4;j++)
for (q.col[i][j]=0,k=0;k<4;k++)
q.col[i][j]+=col[k][j]*m.col[i][k];
return q;
}
_mat4<T> operator * (_vec4<T>&v) // =a[][]*v[]
{
_vec4<T> q;
int i,j;
for (i=0;i<4;i++)
for (q.dat[i]=0,j=0;j<4;j++)
q.dat[i]+=col[i][j]*v.dar[j];
return q;
}
_mat4<T> operator * (T &c) // =a[][]*c
{
_mat4<T> q;
int i,j;
for (i=0;i<4;i++)
for (j=0;j<4;j++)
q.dat[i]=col[i][j]*c;
return q;
}
_mat4<T> operator / (T &c) // =a[][]/c
{
_mat4<T> q;
int i,j;
for (i=0;i<4;i++)
for (j=0;j<4;j++)
q.dat[i]=divide(col[i][j],c);
return q;
}
_mat4<T> operator *=(_mat4<T>&m){ this[0]=this[0]*m; return *this; };
_mat4<T> operator *=(_vec4<T>&v){ this[0]=this[0]*v; return *this; };
_mat4<T> operator *=(const T &c){ this[0]=this[0]*c; return *this; };
_mat4<T> operator /=(const T &c){ this[0]=this[0]/c; return *this; };
// members
void get(T *a)
{
int i,j,k;
for (k=0,i=0;i<4;i++)
for (j=0;j<4;j++,k++)
a[k]=col[i].dat[j];
}
void set(T *a)
{
int i,j,k;
for (k=0,i=0;i<4;i++)
for (j=0;j<4;j++,k++)
col[i].dat[j]=a[k];
}
};
//---------------------------------------------------------------------------
template <class T> _mat4<T> transpose(const _mat4<T> &m)
{
_mat4<T> q;
int i,j;
for (i=0;i<4;i++)
for (j=0;j<4;j++)
q.col[i][j]=m.col[j][i];
return q;
}
//---------------------------------------------------------------------------
template <class T> _mat4<T> inverse(_mat4<T> &m)
{
T p[3];
_mat4<T> q;
T x,y,z;
int i,j;
// transpose rotation
for (i=0;i<3;i++) for (j=0;j<3;j++) q.col[i][j]=m.col[j][i];
// copy projection
for (i=0;i<4;i++) q.col[i][3]=m.col[i][3];
// convert origin: new_pos = - new_rotation_matrix * old_pos
for (i=0;i<3;i++) for (p[i]=0,j=0;j<3;j++) p[i]+=q.col[j][i]*m.col[3][j];
for (i=0;i<3;i++) q.col[3][i]=-p[i];
return q;
}
//---------------------------------------------------------------------------
template <class T> _mat4<T> lrotx(_mat4<T> &m,T ang)
{
T c=cos(ang),s=sin(ang);
_mat4<T> r=mat4(
1, 0, 0, 0,
0, c, s, 0,
0,-s, c, 0,
0, 0, 0, 1);
r=m*r; return r;
};
//---------------------------------------------------------------------------
template <class T> _mat4<T> lroty(_mat4<T> &m,T ang)
{
T c=cos(ang),s=sin(ang);
_mat4<T> r=mat4(
c, 0,-s, 0,
0, 1, 0, 0,
s, 0, c, 0,
0, 0, 0, 1);
r=m*r; return r;
};
//---------------------------------------------------------------------------
template <class T> _mat4<T> lrotz(_mat4<T> &m,T ang)
{
T c=cos(ang),s=sin(ang);
_mat4<T> r=mat4(
c, s, 0, 0,
-s, c, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1);
r=m*r; return r;
};
//---------------------------------------------------------------------------
template <class T> _mat4<T> rotate(_mat4<T> &m,T ang,_vec3<T> p0,_vec3<T> dp)
{
int i;
T c=cos(ang),s=sin(ang);
_vec3<T> x,y,z;
_mat4<T> a,_a,r=mat4(
1, 0, 0, 0,
0, c, s, 0,
0,-s, c, 0,
0, 0, 0, 1);
// basis vectors
x=normalize(dp); // axis of rotation
y=_vec3<T>(1,0,0); // any vector non parallel to x
if (fabs(dot(x,y))>0.75) y=_vec3<T>(0,1,0);
z=cross(x,y); // z is perpendicular to x,y
y=cross(z,x); // y is perpendicular to x,z
y=normalize(y);
z=normalize(z);
// feed the matrix
for (i=0;i<3;i++)
{
a[0][i]= x[i];
a[1][i]= y[i];
a[2][i]= z[i];
a[3][i]=p0[i];
a[i][3]=0;
} a[3][3]=1;
_a=inverse(a);
r=m*a*r*_a;
return r;
};
//---------------------------------------------------------------------------
template <class T> _mat4<T> grotx(_mat4<T> &m,T ang){ return inverse(lrotx(inverse(m),ang)); };
template <class T> _mat4<T> groty(_mat4<T> &m,T ang){ return inverse(lroty(inverse(m),ang)); };
template <class T> _mat4<T> grotz(_mat4<T> &m,T ang){ return inverse(lrotz(inverse(m),ang)); };
//---------------------------------------------------------------------------
typedef _mat4<float > mat4;
typedef _mat4<double> dmat4;
typedef _mat4<bool > bmat4;
typedef _mat4<int > imat4;
typedef _mat4<DWORD > umat4;
//---------------------------------------------------------------------------
mat4 GLSL_math_test4x4;
//---------------------------------------------------------------------------
In order to understand it or write your own I recommend to see:
And lastly I also used mine dynamic list template so:
List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items
Now the result in OBB
its just box described by its center p0 and half vectors u,v,w. So to obtain OBB of pointcloud PCL just compute:
OBB obb;
pointcloud PCL;
PCL.reset();
PCL.add(...); // here feed points into PCL
obb.compute(PCL);
and that is all.
这篇关于计算 3D 网格边界框对角线长度的简单方法是什么?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
