计算 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)。但是我不知道如何利用这些获得的 min/max 来计算边界框的对角线长度。
有什么帮助吗？
为简单起见，让我们考虑 n 个 3D 点（点云）列表作为输入（而不是网格），这对于多边形网格来说已经足够了。
网格的 "diagonal" 只是网格中 2 个最远点之间的线。这很容易通过简单的 O(n^2) 蛮力搜索来计算（2 嵌套 for 循环记住最远的点）。还有一些更快的方法可以利用点的排序。这里是蛮力示例：
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 类 实施的更多信息，请阅读下面的 link 和 OBB 的源代码。
然而，从评论中我感觉你需要 3D OBB（定向边界框），而不仅仅是对角线。你现在拥有的只是 AABB （轴对齐边界框），它不会给你网格对角线（除非它的幸运方向匹配 AABB 对角线）。
注意 AABB 和 OBB 对角线与网格对角线不同！！！
有很多方法可以使用特征向量、凸包等从蛮力 (~O(n^6)) 到更快地计算 OBB...
我设法将我的 移植到 3D 中。
思路是一样的。在 "all" (m) 中存储最大距离可能 directions/angles（覆盖整个球体而不是 2D 中的圆）将数据从 n 减少到 m。然后只需在计算数据中搜索最小边界体积（而不是二维区域）。
我使用 进行测试并作为起点。
算法：
计算枢轴点p0
它必须在 OBB 的内部点。通常 AABB 中心或平均点就足够了。
计算每个可能方向的距离
有无数种可能的方向，因此我们需要将其限制为 m。 m 越大计算越慢但更准确。为了快速存储和获取这些值，我使用了 cube_map.
它是覆盖单位立方体表面（6 x 正方形边）的 2D 纹理，它由方向向量而不是纹理坐标寻址。
我实现了 2 个函数，它们在纹理数据（存储为一维数组）中的 index 和向量中的 direction 之间进行转换。有关详细信息，请参阅示例中的 cube_map...
点 p 从 p0 在某个方向 dir 的距离 d 计算如下：
d = dot( p-p0 , dir )
因此生成 m 个可能的方向，并针对源点列表中所有点的每个计算距离，并记住最大的一个，然后将其存储到 cube_map 以供后者使用。这是O(m*n)
此处为一帧的存储距离示例（cube_map的内容）：
求最小边界体积
简单地生成一些坐标系（覆盖半球）的所有m旋转。你不需要覆盖整个球体，因为另一半只是否定...
现在对于每个计算体积，通过获取沿其 3 个轴的两个方向的距离并计算形成的盒子的体积并记住最小的一个（轴、距离和体积）。由于舍入和非线性问题，cube_map 中可能存在单位化数据，这会导致 volume = 0（如果 cube_map 在开始时被清除为零），因此请忽略此类体积。
在此之后你应该有你的 OBB 近似值。这里预览 OBB 几个旋转位置：
这有点紧张，因为对于这种对称形状，有无限数量的有效 OBBs 并且在不同的旋转中可以在搜索中首先找到不同的。
提高准确性
简单地搜索附近的几个旋转找到 OBB 近似值并记住最小的一个。这可以递归地完成。但是我懒得实施这个，因为 OBB 结果的当前状态对我来说已经足够了。
这里是C++/GL源码（其余的可以在link上面找到）：
//---------------------------------------------------------------------------
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();
}
};
//---------------------------------------------------------------------------
希望我没有忘记复制一些东西...我想让代码尽可能简单，所以它不是很优化，还有很大的改进空间。使用的数学基于 GLSL，因此您可以使用 GLM。我使用自己的库，因为如果需要，vec 可以在上面的 links 中找到（但需要生成为它的 ~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++,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;
//---------------------------------------------------------------------------
为了自己看懂或者自己写推荐看：
最后我还使用了我的动态列表模板：
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 项预分配 space
现在OBB中的结果
它只是由中心 p0 和半向量 u,v,w 描述的盒子。因此，要获得 OBB 的点云 PCL 只需计算：
OBB obb;
pointcloud PCL;
PCL.reset();
PCL.add(...); // here feed points into PCL
obb.compute(PCL);
仅此而已。
我想计算 3D 网格边界框的对角线长度。使用 C++，我迭代顶点并搜索 X 坐标的 (min,max)、Y 坐标的 (min,max) 和 Z 坐标的 (min,max)。但是我不知道如何利用这些获得的 min/max 来计算边界框的对角线长度。 有什么帮助吗？
为简单起见，让我们考虑 n 个 3D 点（点云）列表作为输入（而不是网格），这对于多边形网格来说已经足够了。
网格的 "diagonal" 只是网格中 2 个最远点之间的线。这很容易通过简单的 O(n^2) 蛮力搜索来计算（2 嵌套 for 循环记住最远的点）。还有一些更快的方法可以利用点的排序。这里是蛮力示例：
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 类 实施的更多信息，请阅读下面的 link 和 OBB 的源代码。
然而，从评论中我感觉你需要 3D OBB（定向边界框），而不仅仅是对角线。你现在拥有的只是 AABB （轴对齐边界框），它不会给你网格对角线（除非它的幸运方向匹配 AABB 对角线）。
注意 AABB 和 OBB 对角线与网格对角线不同！！！
有很多方法可以使用特征向量、凸包等从蛮力 (~O(n^6)) 到更快地计算 OBB...
我设法将我的
思路是一样的。在 "all" (m) 中存储最大距离可能 directions/angles（覆盖整个球体而不是 2D 中的圆）将数据从 n 减少到 m。然后只需在计算数据中搜索最小边界体积（而不是二维区域）。
我使用
算法：
计算枢轴点
p0它必须在 OBB 的内部点。通常 AABB 中心或平均点就足够了。
计算每个可能方向的距离
有无数种可能的方向，因此我们需要将其限制为
m。m越大计算越慢但更准确。为了快速存储和获取这些值，我使用了cube_map.它是覆盖单位立方体表面（6 x 正方形边）的 2D 纹理，它由方向向量而不是纹理坐标寻址。
我实现了 2 个函数，它们在纹理数据（存储为一维数组）中的
index和向量中的direction之间进行转换。有关详细信息，请参阅示例中的cube_map...点
p从p0在某个方向dir的距离d计算如下：d = dot( p-p0 , dir )因此生成
m个可能的方向，并针对源点列表中所有点的每个计算距离，并记住最大的一个，然后将其存储到cube_map以供后者使用。这是O(m*n)此处为一帧的存储距离示例（cube_map的内容）：
求最小边界体积
简单地生成一些坐标系（覆盖半球）的所有
m旋转。你不需要覆盖整个球体，因为另一半只是否定...现在对于每个计算体积，通过获取沿其 3 个轴的两个方向的距离并计算形成的盒子的体积并记住最小的一个（轴、距离和体积）。由于舍入和非线性问题，
cube_map中可能存在单位化数据，这会导致volume = 0（如果 cube_map 在开始时被清除为零），因此请忽略此类体积。在此之后你应该有你的 OBB 近似值。这里预览 OBB 几个旋转位置：
这有点紧张，因为对于这种对称形状，有无限数量的有效 OBBs 并且在不同的旋转中可以在搜索中首先找到不同的。
提高准确性
简单地搜索附近的几个旋转找到 OBB 近似值并记住最小的一个。这可以递归地完成。但是我懒得实施这个，因为 OBB 结果的当前状态对我来说已经足够了。
这里是C++/GL源码（其余的可以在link上面找到）：
//---------------------------------------------------------------------------
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();
}
};
//---------------------------------------------------------------------------
希望我没有忘记复制一些东西...我想让代码尽可能简单，所以它不是很优化，还有很大的改进空间。使用的数学基于 GLSL，因此您可以使用 GLM。我使用自己的库，因为如果需要，vec 可以在上面的 links 中找到（但需要生成为它的 ~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++,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;
//---------------------------------------------------------------------------
为了自己看懂或者自己写推荐看：
最后我还使用了我的动态列表模板：
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 项预分配 space
现在OBB中的结果
它只是由中心 p0 和半向量 u,v,w 描述的盒子。因此，要获得 OBB 的点云 PCL 只需计算：
OBB obb;
pointcloud PCL;
PCL.reset();
PCL.add(...); // here feed points into PCL
obb.compute(PCL);
仅此而已。