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 p0
it 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 the m the slower computation but more accurate. In order to store and obtain these values fast I used cube_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 index in texture data (stored as 1D array) and direction vector. For more info see cube_map in the example...
The distance d of point p from p0 in some direction dir is computed like this:
d = dot( p-p0 , dir )
so generate m possible directions, and for each compute distance for all points in your source list of point and remember the biggest one which is then stored to cube_map for latter. This is O(m*n)
Here example of stored distances for one frame (content of cube_map):
find minimal bounding volume
Simply generate all m rotations 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_map which results in volume = 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 li
