-
计算 OBB
所以要么使用 PCA 或这个
获取3DOBB。链接中的代码必须移植到3D，但原理是一样的。这里我的更多 advanced 3D OBB approximation 在 cube_map 上使用递归搜索（这里的代码和方法不如它）。
-
初步猜测
所以 OBB 会给你定向边界框。它最大的一侧将平行于气缸的旋转轴。因此，让我们从描述这个 OBB 的圆柱体开始。所以中心轴将是 OBB 的中心并且平行于它的最大边。 （如果你没有最大的一面，那么你需要检查所有 3 个组合）。直径将是其余边中较大的一个。
-
装缸
现在只需尝试“所有”偏移和半径组合（也可能是高度），将所有点包围在初始猜测附近，并记住最好的一个（根据您想要的规格）。您可以为此使用任何优化方法，但我最喜欢的是：
结果的有效性取决于拟合过程。但是不要太疯狂地使用嵌套拟合，因为复杂性也会变得疯狂。
[Edit1] C++ 中的 3D OBB
我很好奇，今天有时间，所以我对 3D OBB 进行了编码，类似于上面链接的 2D 示例。看起来像它的工作。这是预览：
我使用 2 个集群来验证方向...这里是代码（以简单 C++ 类的形式）：
//---------------------------------------------------------------------------
class OBB3D
{
public:
double p0[3],u[3],v[3],w[3]; // origin,3 axises sorted by size asc
double V,l[3]; // volume, and { |u|,|v|,|w| }
double p[8*3]; // corners
OBB3D() {}
OBB3D(OBB3D& a) { *this=a; }
~OBB3D() {}
OBB3D* operator = (const OBB3D *a) { *this=*a; return this; }
//OBB3D* operator = (const OBB3D &a) { ...copy... return this; }
void compute(double *pnt,int num) // pnt[num] holds num/3 points
{
OBB3D o; // temp OBB values
int i,j;
double a,_a,a0,a1,da,ca,sa; int ea; // search angles
double b,_b,b0,b1,db,cb,sb; int eb;
double c,_c,c0,c1,dc,cc,sc; int ec;
double u0[3],v0[3],pmin[3],pmax[3],q,*qq;
const double deg=M_PI/180.0;
p0[0]=0.0; u[0]=1.0; v[0]=0.0; w[0]=0.0; l[0]=0.0; V=-1.0;
p0[1]=0.0; u[1]=0.0; v[1]=1.0; w[1]=0.0; l[1]=0.0;
p0[2]=0.0; u[2]=0.0; v[2]=0.0; w[2]=1.0; l[2]=0.0;
if (num<3) { V=0.0; return; }
a0=0; a1=360.0*deg; da=10.0*deg; _a=a0;
b0=0; b1= 90.0*deg; db=10.0*deg; _b=b0;
c0=0; c1= 90.0*deg; dc=10.0*deg; _c=c0;
// recursively increase precision
for (j=0;j<5;j++)
{
// try all 3D directions with some step da,db
for (ea=1,a=a0;ea;a+=da){ if (a>=a1) { a=a1; ea=0; } ca=cos(a); sa=sin(a);
for (eb=1,b=b0;eb;b+=db){ if (b>=b1) { b=b1; eb=0; } cb=cos(b); sb=sin(b);
// spherical to cartesian direction
o.w[0]=cb*ca;
o.w[1]=cb*sa;
o.w[2]=sb;
// init u,v from cross product
vector_ld(u0,1.0,0.0,0.0);
if (fabs(vector_mul(u0,o.w))>0.75) // |dot(u,w)>0.75| measn near (anti)parallel
vector_ld(u0,0.0,1.0,0.0);
vector_mul(v0,o.w,u0); // v0 = cross(w,u0)
vector_mul(u0,v0,o.w); // u0 = cross(v0,w)
vector_one(u0,u0); // u0/=|u0|
vector_one(v0,v0); // v0/=|v0|
// try all rotations within u0,v0 plane
for (ec=1,c=c0;ec;c+=dc){ if (c>=c1) { c=c1; ec=0; } cc=cos(c); sc=sin(c);
for (i=0;i<3;i++)
{
o.u[i]=(u0[i]*cc)-(v0[i]*sc);
o.v[i]=(u0[i]*sc)+(v0[i]*cc);
}
// now u,v,w holds potential obb axises socompute min,max
pmin[0]=pmax[0]=vector_mul(pnt,o.u); // dot(pnt,u);
pmin[1]=pmax[1]=vector_mul(pnt,o.v); // dot(pnt,v);
pmin[2]=pmax[2]=vector_mul(pnt,o.w); // dot(pnt,w);
for (i=0;i<num;i+=3)
{
q=vector_mul(pnt+i,o.u); if (pmin[0]>q) pmin[0]=q; if (pmax[0]<q) pmax[0]=q;
q=vector_mul(pnt+i,o.v); if (pmin[1]>q) pmin[1]=q; if (pmax[1]<q) pmax[1]=q;
q=vector_mul(pnt+i,o.w); if (pmin[2]>q) pmin[2]=q; if (pmax[2]<q) pmax[2]=q;
}
// compute V,l from min,max
for (o.V=1.0,i=0;i<3;i++) { o.l[i]=pmax[i]-pmin[i]; o.V*=o.l[i]; }
// remember best solution u,v,w,V,l and compute p0
if ((V<0.0)||(V>o.V))
{
*this=o; _a=a; _b=b; _c=c;
for (i=0;i<3;i++) p0[i]=(pmin[0]*u[i])+(pmin[1]*v[i])+(pmin[2]*w[i]);
}
}
}}
a0=(_a-0.5*da); a1=a0+da; da*=0.1;
b0=(_b-0.5*db); b1=b0+db; db*=0.1;
c0=(_c-0.5*dc); c1=c0+dc; dc*=0.1;
}
// sort axises
{ i=0; qq=u; } // w is max
if (l[1]>l[i]){ i=1; qq=v; }
if (l[2]>l[i]){ i=2; qq=w; }
for (j=0;j<3;j++) { q=w[j]; w[j]=qq[j]; qq[j]=q; } q=l[2]; l[2]=l[i]; l[i]=q;
{ i=0; qq=u; } // v is 2nd max
if (l[1]>l[i]){ i=1; qq=v; }
for (j=0;j<3;j++) { q=v[j]; v[j]=qq[j]; qq[j]=q; } q=l[1]; l[1]=l[i]; l[i]=q;
// compute corners from p0,u,v,w,l
for (i=0;i<3;i++)
{
j=i;
p[j]=p0[i] ; j+=3;
p[j]=p0[i]+(l[0]*u[i]) ; j+=3;
p[j]=p0[i]+(l[0]*u[i])+(l[1]*v[i]) ; j+=3;
p[j]=p0[i] +(l[1]*v[i]) ; j+=3;
p[j]=p0[i] +(l[2]*w[i]); j+=3;
p[j]=p0[i]+(l[0]*u[i]) +(l[2]*w[i]); j+=3;
p[j]=p0[i]+(l[0]*u[i])+(l[1]*v[i])+(l[2]*w[i]); j+=3;
p[j]=p0[i] +(l[1]*v[i])+(l[2]*w[i]); j+=3;
}
}
void gl_draw()
{
glBegin(GL_LINES);
glVertex3dv(p+ 0); glVertex3dv(p+ 3);
glVertex3dv(p+ 3); glVertex3dv(p+ 6);
glVertex3dv(p+ 6); glVertex3dv(p+ 9);
glVertex3dv(p+ 9); glVertex3dv(p+ 0);
glVertex3dv(p+12); glVertex3dv(p+15);
glVertex3dv(p+15); glVertex3dv(p+18);
glVertex3dv(p+18); glVertex3dv(p+21);
glVertex3dv(p+21); glVertex3dv(p+12);
glVertex3dv(p+ 0); glVertex3dv(p+12);
glVertex3dv(p+ 3); glVertex3dv(p+15);
glVertex3dv(p+ 6); glVertex3dv(p+18);
glVertex3dv(p+ 9); glVertex3dv(p+21);
glEnd();
}
} obb;
//---------------------------------------------------------------------------
您只需使用点云数据调用计算，其中 num 是点数的 3 倍。结果存储为单位基向量 u,v,w 和原点 p0 以及每个轴的大小 l[] 或 OBB p
的 8 个角点
只需尝试“所有”球面位置并在w 轴上进行一些步骤，然后尝试垂直于每个极轴位置的所有u,v 和w 记住最小体积OBB。然后以较小的步长仅递归搜索找到的最佳解决方案附近的位置以提高准确性。
我认为这应该提供良好的起点。如果你实现了最小圆而不是u,v 旋转（循环for (ec=1,c=c0;ec;c+=dc)），那么你可以直接从这个搜索中获得你的圆柱体。
代码尚未优化（某些部分如w 轴检查）可以移动到嵌套for循环的较低层。但我想尽可能保持简单易懂。
[Edit2] C++ 中的 3D OBC
我设法修改了我的 3D OBB，方法是用最小封闭圆替换 U,V 搜索（希望我实现它正确，但看起来像......）找到投影在 @987654342 上的所有点的最小封闭 2D 圆@plane，使其成为平行于W 的定向边界圆柱体。我使用了pdf from your link 的第一种方法（使用平分线）。结果如下：
蓝色是3D OBB，棕色/橙色是找到的3D OBC。代码如下：
class OBC3D // 3D Oriented Bounding Cylinder
{
public:
double p0[3],u[3],v[3],w[3]; // basecenter,3 axises
double V,r,h; // volume, radius height
double p1[3]; // other base center
OBC3D() {}
OBC3D(OBC3D& a) { *this=a; }
~OBC3D() {}
OBC3D* operator = (const OBC3D *a) { *this=*a; return this; }
//OBC3D* operator = (const OBC3D &a) { ...copy... return this; }
void compute(double *pnt,int num) // pnt[num] holds num/3 points
{
OBC3D o; // temp OBB values
int i,j,k,kk,n;
double a,_a,a0,a1,da,ca,sa; int ea; // search angles
double b,_b,b0,b1,db,cb,sb; int eb;
double pmin[3],pmax[3],q,qq,*pnt2,p[3],c0,c1,u0,v0,du,dv,dr;
const double deg=M_PI/180.0;
p0[0]=0.0; u[0]=1.0; v[0]=0.0; w[0]=0.0; V=-1.0;
p0[1]=0.0; u[1]=0.0; v[1]=1.0; w[1]=0.0; r=0.0;
p0[2]=0.0; u[2]=0.0; v[2]=0.0; w[2]=1.0; h=0.0;
if (num<3) { V=0.0; return; }
// prepare buffer for projected points
pnt2=new double[num];
a0=0; a1=360.0*deg; da=10.0*deg; _a=a0;
b0=0; b1= 90.0*deg; db=10.0*deg; _b=b0;
// recursively increase precision
for (k=0;k<5;k++)
{
// try all 3D directions with some step da,db
for (ea=1,a=a0;ea;a+=da){ if (a>=a1) { a=a1; ea=0; } ca=cos(a); sa=sin(a);
for (eb=1,b=b0;eb;b+=db){ if (b>=b1) { b=b1; eb=0; } cb=cos(b); sb=sin(b);
// spherical to cartesian direction
o.w[0]=cb*ca;
o.w[1]=cb*sa;
o.w[2]=sb;
// init u,v from cross product
vector_ld(o.u,1.0,0.0,0.0);
if (fabs(vector_mul(o.u,o.w))>0.75) // |dot(u,w)>0.75| measn near (anti)parallel
vector_ld(o.u,0.0,1.0,0.0);
vector_mul(o.v,o.w,o.u); // v0 = cross(w,u0)
vector_mul(o.u,o.v,o.w); // u0 = cross(v0,w)
vector_one(o.u,o.u); // u0/=|u0|
vector_one(o.v,o.v); // v0/=|v0|
// now u,v,w holds potential obb axises so compute min,max and convert to local coordinates
pmin[0]=pmax[0]=vector_mul(pnt,o.u); // dot(pnt,u);
pmin[1]=pmax[1]=vector_mul(pnt,o.v); // dot(pnt,v);
pmin[2]=pmax[2]=vector_mul(pnt,o.w); // dot(pnt,w);
for (i=0;i<num;i+=3)
{
q=vector_mul(pnt+i,o.u); if (pmin[0]>q) pmin[0]=q; if (pmax[0]<q) pmax[0]=q; pnt2[i+0]=q;
q=vector_mul(pnt+i,o.v); if (pmin[1]>q) pmin[1]=q; if (pmax[1]<q) pmax[1]=q; pnt2[i+1]=q;
q=vector_mul(pnt+i,o.w); if (pmin[2]>q) pmin[2]=q; if (pmax[2]<q) pmax[2]=q; pnt2[i+2]=q;
}
// [compute min enclosing circle]
n=0;
// center (u0,v0) = avg( pnt2 )
for (u0=0.0,v0=0.0,i=0;i<num;i+=3)
{
u0+=pnt2[i+0];
v0+=pnt2[i+1];
} q=3.0/double(num); u0*=q; v0*=q;
// r = max(|pnt2 - (u0,v0)|)
for (o.r=0.0,i=0;i<num;i+=3)
{
c0=pnt2[i+0]-u0;
c1=pnt2[i+1]-v0;
q=(c0*c0)+(c1*c1);
if (o.r<q) o.r=q;
} o.r=sqrt(o.r);
for (kk=0;kk<4;kk++)
{
// update edgepoints count n
qq=o.r*o.r;
for (i=n;i<num;i+=3)
{
c0=pnt2[i+0]-u0;
c1=pnt2[i+1]-v0;
q=fabs((c0*c0)+(c1*c1)-qq);
if (q<1e-10)
{
pnt2[n+0]=pnt2[i+0];
pnt2[n+1]=pnt2[i+1];
pnt2[n+2]=pnt2[i+2]; n+=3;
}
}
// compute bisector (du,dv)
for (du=0.0,dv=0.0,i=0;i<n;i+=3)
{
du+=pnt2[i+0]-u0;
dv+=pnt2[i+1]-v0;
} q=1.0/sqrt((du*du)+(dv*dv)); du*=q; dv*=q;
// try to move center towards edge points as much as possible
for (dr=0.1*o.r,j=0;j<5;)
{
u0+=dr*du;
v0+=dr*dv;
// q = max(|pnt2 - (u0,v0)|)
for (qq=0.0,i=0;i<num;i+=3)
{
c0=pnt2[i+0]-u0;
c1=pnt2[i+1]-v0;
q=(c0*c0)+(c1*c1);
if (qq<q) qq=q;
} qq=sqrt(qq);
// recursively increase precision
if (qq>o.r)
{
u0-=dr*du;
v0-=dr*dv;
dr*=0.1;
j++;
}
else o.r=qq;
}
}
// compute h,V
o.h=pmax[2]-pmin[2];
o.V=M_PI*o.r*o.r*o.h;
// remember best solution u,v,w,V,l and compute p0
if ((V<0.0)||(V>o.V))
{
*this=o; _a=a; _b=b;
for (i=0;i<3;i++) p0[i]=(u0*u[i])+(v0*v[i])+(pmin[2]*w[i]);
}
}}
a0=(_a-0.5*da); a1=a0+da; da*=0.1;
b0=(_b-0.5*db); b1=b0+db; db*=0.1;
}
// compute corners from p0,u,v,w,l
for (i=0;i<3;i++) p1[i]=p0[i]+(h*w[i]);
delete[] pnt2;
}
void gl_draw()
{
int i,j,n=36;
double a,da=2.0*M_PI/double(n),p[3],uu,vv;
glBegin(GL_LINES);
glVertex3dv(p0); glVertex3dv(p1);
glEnd();
glBegin(GL_LINE_LOOP);
for (a=0.0,i=0;i<n;i++,a+=da)
{
uu=r*cos(a);
vv=r*sin(a);
for (j=0;j<3;j++) p[j]=p0[j]+(u[j]*uu)+(v[j]*vv);
glVertex3dv(p);
}
glEnd();
glBegin(GL_LINE_LOOP);
for (a=0.0,i=0;i<n;i++,a+=da)
{
uu=r*cos(a);
vv=r*sin(a);
for (j=0;j<3;j++) p[j]=p1[j]+(u[j]*uu)+(v[j]*vv);
glVertex3dv(p);
}
glEnd();
}
};
//---------------------------------------------------------------------------
用法是一样的...我用这个测试过：
OBB3D obb;
OBC3D obc;
void compute()
{
int i,n=500;
// random pnt cloud
Randomize();
RandSeed=98123456789;
pnt.allocate(3*n); pnt.num=0;
// random U,V,W basis vectors
double u[3],v[3],w[3],x,y,z,a;
for (i=0;i<3;i++) w[i]=Random()-0.5; // random direction
vector_one(w,w); // w/=|w|
vector_ld(u,1.0,0.0,0.0);
if (fabs(vector_mul(u,w))>0.75) // |dot(u,w)>0.75| measn near (anti)parallel
vector_ld(u,0.0,1.0,0.0);
vector_mul(v,w,u); // v = cross(w,u)
vector_mul(u,v,w); // u = cross(v,w)
vector_one(u,u); // u/=|u|
vector_one(v,v); // v/=|v|
// random cylinder point cloud
for (i=0;i<n;i++)
{
a=2.0*M_PI*Random();
x= 0.5+(0.75*(Random()-0.5))*cos(a);
y=-0.3+(0.50*(Random()-0.5))*sin(a);
z= 0.4+(0.90*(Random()-0.5));
pnt.add((x*u[0])+(y*v[0])+(z*w[0]));
pnt.add((x*u[1])+(y*v[1])+(z*w[1]));
pnt.add((x*u[2])+(y*v[2])+(z*w[2]));
}
obb.compute(pnt.dat,pnt.num);
obc.compute(pnt.dat,pnt.num);
}
List<double> pnt 是我的动态数组模板double pnt[]。在这里哪个不重要。
请注意，如果您为 W 方向搜索选择了太大的初始步长 (da,db)，您可能会因为将自身困在局部最小值内而错过正确的解决方案。