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.