我想计算3D网格的边界框的对角线长度。使用C ++迭代顶点,并搜索X坐标的(min,max),Y坐标的(min,max)和Z坐标的(min,max)。但是,我不知道如何利用这些获得的最小值/最大值来计算边界框的对角线长度。有什么帮助吗?
为简单起见,我们考虑将n
3D点(点云)的列表作为输入(而不是网格),足以用于多边形网格。
网格的“对角线”就是网格中两个最远点之间的直线。这可以通过简单的O(n^2)
蛮力搜索(2个嵌套循环来记住最远的点)轻松地计算出来。还有一些更快的方法可以利用点的排序。这里是蛮力的例子:
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;
}
有关line
和pointcloud
类实现的更多信息,请阅读下面的链接和OBB的源代码。
但是从评论中我感觉到您需要3D OBB(定向边界框),而不只是对角线。您现在拥有的只是AABB(与轴对齐的边界框),它不会给您网格对角线(除非它的幸运方向与AABB对角线匹配)。
当心AABB和OBB对角线与网格对角线不同!!
有很多方法可以使用特征向量,凸包等从蛮力(〜O(n^6)
)更快地计算OBB ...
我设法将2D OBB approximation移植到3D中。
想法是一样的。将最大距离存储在“所有”(m
)可能的方向/角度(在2D中覆盖整个球面而不是圆形),以将数据从n
减少到m
。然后只搜索计算出的数据以获取最小边界体积(而不是2D中的面积)。
我将Cone to box collision用于测试并作为起点。
算法:
计算枢轴点p0
它必须在OBB的内部。通常,AABB的中心或平均点就足够了。
计算每个可能方向的距离
存在无限可能的方向,因此我们需要将其限制为m
。 m
越大,计算速度越慢,但精度更高。为了快速存储和获取这些值,我使用了cube_map
。
它是一个二维纹理,覆盖单位立方体的表面(6个正方形),并通过方向矢量而不是纹理坐标进行寻址。
我实现了2个函数,用于在纹理数据(存储为1D数组)中的index
和direction
向量之间进行转换。有关更多信息,请参见示例中的cube_map
...
d
与p
在某个方向p0
上的距离dir
的计算如下:
d = dot( p-p0 , dir )
因此生成m
个可能的方向,并针对点源列表中所有点的每个计算距离,记住最大的一个,然后将其存储到cube_map
。这是O(
m * n )
这里存储一帧距离的示例(cube_map的内容):
寻找最小边界体积
仅生成某个坐标系(覆盖半球)的所有m
旋转。您不需要覆盖整个领域,因为另一半只是否定...
现在,通过获取两个方向上沿其3个轴的距离并计算形成的框的体积,并记住最小的框(轴,距离和体积)来计算每个体积。由于舍入和非线性问题,在cube_map
中可能存在统一的数据,从而导致volume = 0
(如果cube_map在开始时清零),因此请忽略此类体积。
在此之后,您应该使用OBB代替。这是OBB的几个旋转位置的预览:
这有点跳跃,因为对于这种对称形状,有无限数量的有效OBB,并且在不同的旋转中,可以首先在搜索中找到不同的一个。]]
提高准确性
简单地搜索附近发现的几个旋转<< OBB
并记住最小的旋转。这可以递归完成。但是,我实在太懒惰了,因为OBB结果的当前状态对我来说足够了。此处为C ++ / GL源(其余内容可在上面的链接中找到:
//---------------------------------------------------------------------------
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();
}
};
//---------------------------------------------------------------------------
希望我没有忘记复制某些东西……我想使代码尽可能地简单,因此它不是很优化,还有很多改进的余地。使用的数学基于,因此您可以使用GLM。我使用了自己的库,以便在需要时可以在上面的链接中找到GLSL
vec
(但需要将其作为〜220KByte的代码生成),但是它与GLSL和GLM完全匹配,因此您可以使用cn。但是mat4
使用的某些功能在GLM中不存在,因此以防万一: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++)
a[k]=col[i].dat[j];
}
};
//---------------------------------------------------------------------------
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(const _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]*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 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 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 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> 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;
//---------------------------------------------------------------------------
为了理解它或自己写,我建议看看:
List<double> xxx;
与double xxx[];
相同xxx.add(5);
将5
添加到列表的末尾xxx[7]
访问数组元素(安全)xxx.dat[7]
访问数组元素(不安全但快速的直接访问)xxx.num
是数组的实际使用大小xxx.reset()
清除数组并设置xxx.num=0
xxx.allocate(100)
为100
个项目预分配空间
现在将结果放入OBB
仅由其中心p0
和一半向量u,v,w
描述的方框。因此,要获取点云PCL
的OBB
,只需计算:OBB obb;
pointcloud PCL;
PCL.reset();
PCL.add(...); // here feed points into PCL
obb.compute(PCL);
仅此而已。