计算 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(){整数 i, j;线 l,ll;l=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0));//空行for (i=0;i有关 line 和 pointcloud 类实现的更多信息，请阅读以下链接和 OBB 的源代码.
但是从评论中我感觉你需要3D OBB(定向边界框)而不是对角线.您现在拥有的只是 AABB(轴对齐边界框)，它不会为您提供网格对角线(除非它的幸运方向与 AABB 对角线匹配).
注意 AABB 和 OBB 对角线与网格对角线不同！！！
有许多方法可以使用特征向量、凸包等从蛮力 (~O(n^6)) 更快地计算 OBB...
我设法将我的
找到最小包围体
只需生成某个坐标系(覆盖半个球体)的所有 m 旋转.你不需要覆盖整个领域，因为另一半只是否定......
现在通过获取沿其 3 个轴在两个方向上的距离并计算形成的盒子的体积并记住最小的一个(轴、距离和体积)来计算每个体积.由于舍入和非线性问题，可能会在 cube_map 中包含未初始化的数据，这会导致 volume = 0(如果 cube_map 在开始时被清除为零)，因此请忽略这些只是体积.
在此之后，您应该有您的 OBB 近似值.这里是几个旋转位置的 OBB 预览:
有点跳跃，因为对于这种对称的形状，有无数个有效的OBB，并且在不同的旋转中，可以在搜索中首先找到不同的.
提高准确性
简单地在附近搜索几圈找到 OBB 近似值并记住最小的那个.这可以递归地完成.但是我懒得实现这一点，因为 OBB 结果的当前状态对我来说已经足够了.
这里是 C++/GL 源码(其余的可以在上面的链接中找到):
//--------------------------------------------------------------------------类点云{上市://配置文件列出<vec3>点；点云(){}点云(点云& a){ *this=a;}〜点云(){}pointcloud* operator = (const pointcloud *a) { *this=*a;返回这个；}//pointcloud* operator = (const pointcloud &a) { ...copy... return this;}无效重置(){ pnt.num = 0;}无效添加(vec3 p){ pnt.add(p)；}无效添加(点 p){ pnt.add(p.p0);}无效计算(){}；无效绘制(){glBegin(GL_POINTS);for (int i=0;i=sz)) 返回目录；x=ix%n;ix/=n；x/=fn2；X - ;y=ix%n;ix/=n；y/=fn2；是——；如果(ix==0){ dir.y=x;dir.z=y;dir.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;目录.z=-1.0;}如果(ix==5){ dir.x=x;目录.y=y;目录.z=+1.0；}返回规范化(目录)；}int dir2ix(vec3 目录){int ix=0,x=0,y=0;浮动 a=0.0,b;b=fabs(dir.x);如果 (a=sz)) ix=0;返回九；}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未旋转坐标系(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;aVV)&&(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;}无效绘制(){常量 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，};常量 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 中不存在的函数，以防万一:
模板<class T>类_mat4{上市:_vec4<T>科尔[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* 运算符 = (const _mat4 &a) { for (int i=0;i<4;i++) col[i]=a.col[i];返回这个；}//=a[][]_vec4<T>&运算符 [](const int i){ return col[i];}//a[i]_mat4<T>运算符 * (_mat4&m)//=a[][]*m[][]{_mat4<T>q;整数 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<T>运算符 * (_vec4&v)//=a[][]*v[]{_vec4<T>q;整数 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>运算符 * (T &c)//=a[][]*c{_mat4<T>q;整数 i, j;对于 (i=0;i<4;i++)对于 (j=0;j<4;j++)q.dat[i]=col[i][j]*c;返回q；}_mat4<T>运算符/(T &c)//=a[][]/c{_mat4<T>q;整数 i, j;对于 (i=0;i<4;i++)对于 (j=0;j<4;j++)q.dat[i]=除法(col[i][j],c);返回q；}_mat4<T>运算符 *=(_mat4&m){ this[0]=this[0]*m;返回*这个；};_mat4<T>运算符 *=(_vec4&m){_mat4<T>q;整数 i, j;对于 (i=0;i<4;i++)对于 (j=0;j<4;j++)q.col[i][j]=m.col[j][i];返回q；}//--------------------------------------------------------------------------模板<T类>_mat4<T>逆(_mat4&m){T p[3];_mat4<T>q;T x,y,z;整数 i, j;//转置旋转对于 (i=0;i<3;i++) 对于 (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_pos对于 (i=0;i<3;i++) 对于 (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；}//--------------------------------------------------------------------------模板<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;返回 r;};//--------------------------------------------------------------------------模板<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;返回 r;};//--------------------------------------------------------------------------模板<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;返回 r;};//--------------------------------------------------------------------------模板<T类>_mat4<T>旋转(_mat4p0,_vec3dp){诠释我;T c=cos(ang),s=sin(ang);_vec3;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);//基向量x=标准化(dp)；//旋转轴y=_vec3(1,0,0);//任何与 x 不平行的向量if (fabs(dot(x,y))>0.75) y=_vec3<T>(0,1,0);z=交叉(x，y)；//z 垂直于 x,yy=交叉(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;};//--------------------------------------------------------------------------模板<T类>_mat4<T>grotx(_mat4<T> &m,T ang){ return inverse(lrotx(inverse(m),ang));};模板<T类>_mat4<T>groty(_mat4<T> &m,T ang){ return inverse(lroty(inverse(m),ang));};模板<T类>_mat4<T>grotz(_mat4<T> &m,T ang){ return inverse(lrotz(inverse(m),ang));};//--------------------------------------------------------------------------typedef _mat4垫4;typedef _mat4<double>dmat4;typedef _mat4bmat4;typedef _mat4易马4;typedef _mat4乌马特4；//--------------------------------------------------------------------------mat4 GLSL_math_test4x4；//-------------------------------------------------------------------------- 为了理解或者自己写我推荐看看:
最后我还使用了我的动态列表模板:
列表<double>xxx; 等同于 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屋！
