最小外接圆柱

Question

是否有一种算法可以找到一个能够用最小半径包围3D点云的圆柱体？我知道2D情况下可以解决最小包围圆问题（例如此线程 Python中的最小包围圆，代码错误），但是在3D下是否有任何可行的方法？

---

编辑1: OBB。以下是一个弧形点云的示例。该工具https://www.nayuki.io/page/smallest-enclosing-circle 找到了最小包围圆。

圆由三个点定义，其中两个点几乎位于直径上，因此很容易估计中心轴的位置。 "打包"点将得到一个明显偏离真实中心的盒子中心。

我得出结论： OBB方法不是通用的。

---

编辑2: PCA. 下面是一个关于紧密点云与带有离群值的点云的PCA分析示例。对于紧密点云，PCA能够令人满意地预测圆柱方向。但是如果与主点云相比，存在少量离群点，则PCA基本上会忽略它们，得出的向量距离一个包围圆柱的真实轴线非常远。在下面的示例中，一个包围圆柱的真实几何轴线以黑色显示。

我认为PCA方法不具有普适性。

---

编辑3: OBB vs. PCA和OLS. 一个重要的区别 - OBB仅依赖于几何形状，而PCA和OLS则依赖于整个点数，包括不影响形状的中间集合中的那些点。为了使它们更有效，可以包括数据准备步骤。首先找到凸壳。其次，排除所有内部点。然后，沿着凸壳的点可以不均匀分布。我建议删除所有这些点，只留下多边形凸包体，并用网格覆盖它，其中节点将成为新的点。对这些新点云应用PCA或OLS应该可以更准确地估计圆柱轴。

如果OBB提供了一个与外包圆柱轴尽可能平行的轴，那么所有这些都是不必要的。

---

编辑4：已发布的方法。

@meowgoesthedog：Michel Petitjean的论文（“关于最小外接圆柱问题的代数解法”）可能有所帮助，但我不够资格将其转换为可工作的程序。作者自己做到了（模块CYL在此处http://petitjeanmichel.free.fr/itoweb.petitjean.freeware.html）。但在论文的结论中，他说：“而且这个名为CYL的现有软件，可在http://petitjeanmichel.free.fr/itoweb.petitjean.freeware.html免费下载，既不声称提供最佳方法的实现，也不声称比其他圆柱计算软件更好。”论文中的其他短语也给人留下印象，即它是一种实验性方法，没有经过彻底验证。无论如何，我将尝试使用它。

@Ripi2：Timothy M. Chan的这篇论文对我来说也有点复杂。我不是数学方面的专家，无法将其转换为工具。

@Helium_1s2：可能这是一个好建议，但与上述两篇论文相比，它要少得多。而且没有经过验证。

编辑5：为用户1717828回复。两个最远的点与圆柱轴线。

反例 - 一个立方体的8个点，适合于圆柱体内。最大的距离在两点之间的绿色对角线。显然不与圆柱轴平行。

"中点"方法由Ripi2提出：仅适用于二维情况。在三维情况下，圆柱轴可能不与任何两点之间的单线段相交。

- JRP

---

这个问题让人想起了PCA。 - collapsar

那么你需要查看这篇论文。https://www.google.com.ar/url?sa=t&source=web&rct=j&url=https://www.geometrictools.com/Documentation/CylinderFitting.pdf&ved=2ahUKEwiV98C30K3cAhUbfisKHYIIAjYQFjABegQIAhAB&usg=AOvVaw0FVGNisCPiA4aPceCJ7F17 - Mauricio Cele Lopez Belon

@JRP，我刚刚看到你编辑了你的问题+1。3D OBB和PCA应该导致相同的结果，但你是对的，这不是你的解决方案，因为最佳候选者可能会稍微偏离找到的轴、中心和大小，但不会太多，这就是为什么之后需要使用拟合。然而，拟合可能会有点简化，因为你可以直接从数据集计算出圆柱的高度，只留下2D轴和半径要拟合的偏移量。半径也可以从数据集计算出来，所以只剩下底部中心的x，y要拟合，这应该不难。你有样本数据集吗？ - Spektre

@Spektre: 我怀疑 OBB 和 PCA 不会给出相同的结果 - 我会将其添加到问题中。任何随机数据集都将适合此时。 - JRP

@JRP 我在我的答案中添加了一个C++示例，名为Edit2，其中包含了3D OBC。看起来它正在工作，但如何证明它是个问题... - Spektre

Answer 1

英译中：

1. 计算OBB

   可以使用PCA或者以下方法：

   • 如何计算多条曲线的OBB？

   来获取3D OBB。链接中的代码需要转换成3D，但原理是相同的。这里有一个更高级的3D OBB逼近，使用了立方体映射的递归搜索（这里的代码和方法不如它）。

2. 初始猜测

   因此，OBB将给出定向边界框。其最大的一面将与您的圆柱体的旋转轴平行。因此，让我们从圆柱体外围绕绕着这个OBB开始。因此，中心轴将是OBB的中心，并且与其最大的一面平行。（如果没有最大的一面，则需要检查所有3种组合）。直径将是剩余面中的较大者。

3. 拟合圆柱体

   现在只需尝试囊括初始猜测附近所有点的“所有”偏移量和半径（也可能是高度），并记住最佳的一个（根据您想要的规格）。可以使用任何优化方法进行此操作，但我最喜欢的是这个：

   • 逼近搜索如何工作

结果的有效性取决于拟合过程。但是不要在嵌套拟合中变得过于复杂。

[编辑1] C++中的3D OBB

我很好奇，今天有些时间，所以我编码了类似于上面链接的2D示例的3D OBB。看起来它正在工作。这是预览：

我用了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的8个角点p。

这个算法的实现非常简单，通过尝试所有球形位置（w轴有一定的步长），然后尝试所有与w垂直的u,v极坐标位置，并记住最小体积的OBB。然后递归搜索附近的最佳解决方案，以提高精度。

我认为这应该是一个不错的起点。如果您在循环for (ec=1,c=c0;ec;c+=dc)中实现最小圆而不是u,v旋转，则可以直接从此搜索中获得您的圆柱体。

代码还没有优化（例如w轴检查等部分可以移动到嵌套的for循环的较低层），但我想尽可能保持它简单易懂。

我成功地修改了我的3D OBB，通过用最小外接圆替换U,V搜索（希望我实现得正确，但看起来是这样的...），找到所有投影在UV平面上的点的最小外接2D圆，使其成为与W平行的定向包围圆柱。我使用了你提供的链接中的第一种方法（使用分界线）。以下是结果：

蓝色部分是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[]，但这并不重要。

请注意，如果您选择了过大的初始步长（da,db）进行 W 方向搜索，可能会陷入局部最小值而错过正确的解。

- Spektre

Answer 2

概念性答案

1. 找到距离最远的两个点h。它们在圆柱体的面上，连接它们的线将与圆柱体的轴平行。

2. 投影所有点到垂直于该轴的平面上。

3. 找到平面上距离最远的两个点d。它们定义了直径为d的圆。

4. 包含所有点的最小可能体积的圆柱体为

.

* 这假设只有一对定义圆柱体轴的距离最远的两个点。如果有可能有两对点共享最大值，请针对每一对重复步骤2-4，并选择直径最小的圆柱体。

Python答案

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook

from numpy.linalg import norm
from scipy.spatial.distance import pdist, squareform

如果您还没有生成点，请执行以下操作：

np.random.seed(0)
N = 30
M = np.random.randint(-3,3,(N,3))
print(M)

[[ 1  2 -3]
 [ 0  0  0]
 [-2  0  2]
 [-1  1 -3]
 [-3  1 -1]
 ...
 [ 1 -3  1]
 [ 0 -1  2]]

计算每对点之间的距离，并选择距离最大的一对。

max_dist_pair = list(pd.DataFrame(squareform(pdist(M))).stack().idxmax())
p1 = M[max_dist_pair[0]]
p2 = M[max_dist_pair[1]]
print(f"Points defining cylinder faces: {p1}, {p2}")
print(f"Length of cylinder: {norm(p1-p2)}")

Points defining cylinder faces: [-1 -3 -3], [1 2 2]
Length of cylinder: 7.3484692283495345

将所有点都用蓝色图示，最大间隔的点用红色标出。

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(*M.T, c = ['red' if i in max_dist_pair else 'blue' for i in range(N)])

ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")

plt.show()

这里是同一图形的旋转视角，我们现在看向两个红点之间的轴。

以上视图相当于将点投影到垂直于圆柱体轴的平面上。在该平面上找出最小的包含所有点的圆。我们通过计算每个点到轴的位移，并找出两点之间的最大距离来实现此目标。

perp_disp = (np.cross(p2-p1, M-p1))/norm(p2-p1) # Get perpendicular displacement vectors.
print(perp_disp)

[[-3.40206909  1.36082763  0.        ]
 [ 0.         -0.13608276  0.13608276]
 [ 1.36082763 -2.04124145  1.4969104 ]
 [-2.72165527  0.          1.08866211]
 [-1.36082763 -1.90515869  2.44948974]
 [ 0.68041382 -0.95257934  0.68041382]
 [ 2.72165527  0.68041382 -1.76907593]
 ...
 [ 0.          0.27216553 -0.27216553]
 [ 0.         -0.40824829  0.40824829]
 [ 2.72165527  0.27216553 -1.36082763]
 [ 2.04124145 -0.68041382 -0.13608276]]

通过使用与上面相同的“pdist”技巧，可以找到最大距离。

max_perp_disp_pair = list(pd.DataFrame(squareform(pdist(perp_disp))).stack().idxmax())
perp_p1 = M[max_perp_disp_pair[0]]
perp_p2 = M[max_perp_disp_pair[1]]
print(perp_p1, perp_p2)

[ 1  2 -3] [-3 -2  1]

最后，我们得到了圆柱的直径。

print(norm(perp_p1 - perp_p2))
6.92820323028

可以包含这些点的圆柱体的最小体积为

注释

• 使用Numpy的两两距离函数pdist找到点之间的最大距离。然后使用squareform将其格式化为Pandas DataFrame，以便可以使用idxmax轻松找到两个点的索引。可能有一种更好的方法来避免使用Pandas。

• 如果np.cross让您感到困惑，那么它仅仅是为了找到点与直线之间的最小距离。如果您对此感兴趣，我可以提供更多详细信息，但是如果您画出两条线的叉积，您会得到一个平行四边形，其中非平行边由这两条线给出。该平行四边形与长度等于其中一条线段的矩形面积相同，宽度等于点到该线的距离。

- user1717828

---

请查看顶部的问题 - EDIT5。我需要粘贴一张图片，但无法在评论中进行。 - JRP

Answer 3

首先找到点云的定向边界框（OBB）。有几种算法可以实现这个目标。以下算法可能是最优的：

https://pdfs.semanticscholar.org/a76f/7da5f8bae7b1fb4e85a65bd3812920c6d142.pdf

现在，可以通过将OBB沿着其最长轴旋转来轻松地找到非最优定向圆柱体，同样，可以找到由OBB包围的圆柱体，其轴与另一个相同，但半径为OBB面法线最短边的一半。

我的猜想是最优圆柱的半径介于这两个圆柱之间。

如果计算所有点到外圆柱的最小距离并调整其半径使该距离等于零，则可以轻松找到最佳圆柱。

这种方法可能有效，但在计算上不太优化，因为您必须计算从所有点到圆柱的距离。也许内部圆柱可以用来剔除所有位于其中的点。我没有详细阐述那个想法。

更新：

似乎这个问题不清楚什么是“最小的”，实际上需要超出“最小”的东西，而且不太清晰。 “最小的”圆柱体包含一组点，应该最小化圆柱体内的空间（至少我理解为最小）。但OP还强制要求最小的圆柱体应该适合输入数据的形状。这意味着，如果输入数据采样了半个圆柱体（通过其最长的一侧切割），答案应该是最符合那一半形状的圆柱体。无论那个圆柱体是否比其他包围数据的圆柱体有更多的空余空间。

这两个要求是相互矛盾的。因为最小的圆柱体可能无法适应数据的曲线形状，而最适合数据曲线形状的圆柱体可能不是最小的圆柱体。

我的答案（和其他答案）基于OBB，回答了关于“最小”的圆柱体包含数据并最小化圆柱体内空间的问题。

另一种将圆柱体与数据形状匹配的情况也可以使用优化方法来回答。但是没有通用的答案。 “最佳”的圆柱体取决于应用程序的需求，并且必须使用至少两种不同的策略根据数据进行计算。

- Mauricio Cele Lopez Belon

---

问题的第一句指出需要“最小半径”。标题不太精确，有点类似于“最小圆”问题的参考。我理解了您关于OBB和优化圆柱体积的答案。然而，OBB中最长边的方向可能近似于圆柱轴的位置。如果是这样的话，就可以简化任务。 - JRP

Answer 4

寻找确切解决方案似乎是一个非常困难的问题。通过对轴向进行假设，并旋转云（仅保留凸包的顶点）并将点投影到XY平面上，您确实可以转向最小包围圆问题。

但是您必须对可能的方向进行多个假设。对于确切解决方案，您应该尝试所有方向，以便围绕不同的顶点对或三元组定义包围圆。这在旋转角度空间中定义了复杂的区域，并且对于每个区域都有一个实现最优化的点。这涉及高度非线性的最小化问题和非线性约束，并且似乎难以处理。

在这个阶段，我所能推荐的是一个近似解决方案，即尝试固定数量的预定义方向并解决相应的圆拟合问题。由于解决方案本身就是近似的，因此近似的圆拟合也可以胜任。

- user1196549

---

我同意你的建议，寻找凸包并删除所有内部点。让我感到好奇的是，用OBB轴逼近圆柱轴的精度如何。 - JRP

面向边界框（Oriented Bounding Box）- 请参考Spektre的最佳答案。 - JRP

OBB 的最长边的方向应该与外接圆柱体的轴线非常接近。 - JRP

有一种算法可以在3D点云周围找到OBB - 请参考Spektre的答案和图像（现在在上面）。如果你想象一下在同样的点云周围画一个最小的包容圆柱体，它看起来像这个圆柱体的轴线将接近于与OBB的最长边平行。因此，OBB的最长边的方向可以用于近似圆柱体的轴线。 - JRP

构建OBB不需要对初始方向做任何假设。以下是一些链接：https://en.wikipedia.org/wiki/Minimum_bounding_box_algorithms，http://graphics.stanford.edu/courses/cs468-06-winter/Slides/an_bounding_box_fall.pdf，https://pdfs.semanticscholar.org/a76f/7da5f8bae7b1fb4e85a65bd3812920c6d142.pdf，http://graphics.stanford.edu/courses/cs468-06-winter/Papers/har-peled-bbox.pdf - JRP

Answer 5

暴力算法

让我们从最简单的情况开始：3个点。

最小的圆柱体表面有三个点。最小的半径意味着两个点在截面的直径上，即使该截面不垂直于圆柱体的轴线。第三个点也是同样的情况。
因此，轴线穿过直径的中心，这是由两点定义的线段的中点。
因此，轴由两个中点定义。

我们有三个中点和三个可能的轴:

通过找到最小ri，选择最佳解决方案（最小半径）。

请注意，每个轴ai都与一些Pi,Pj线段平行。

通用情况，n个点

假设我们已经找到了解决方案。至少有三个点位于表面上，类似于3点情况。如果我们找到了这个三元组，那么我们可以应用中间点方法。 因此，解决方案的轴是通过两个中间点的某条线。在这个暴力算法中，我们测试它们所有。对于每个轴，计算所有垂直于该轴的点的距离，并获得最大值di，因此我们获得每个轴的封闭半径。

解决方案是最小di的轴。半径是这个最小的di。

这个算法的顺序是什么？ 对于n个点，我们有C(n,2)= n(n-1)/2个中间点。和n(n-1)(n(n-1)-2)/8个轴。对于每个轴，我们测试n个半径。所以我们得到O(n^5)。

改进

• 首先构建凸包并丢弃所有内部点
• 我有一种感觉，解决方案的轴由最长的Mi,Mj线段定义。

---

编辑

正如@JRP所示，这种方法不能找到棱柱的最优解。将尝试使用三角形中心（而不是线段中点）进行计算也不起作用，想象一个有8个点的立方体顶点的点云。

这可能是一个很好的近似，但不是最优解。

- Ripi2

---

在3D情况下寻找中点可能会存在困难-请参见主问题中的EDIT5，三棱柱情况。 - JRP

@JRP 嗯，看起来我的方法比你画的三角形中心情况给出了更好的解决方案（最小半径）。 - Ripi2

@JRP。你对我的解决方案是正确的。也许使用三角形中心而不是线段中点可能会起作用。 - Ripi2

在这篇论文中（https://arxiv.org/pdf/1008.5259.pdf），他们从4点云（四面体）和5点云（三角双锥体）开始分析3D情况，并成功地通过解析方法找到了轴。 - JRP

但我无法理解他们如何在大量点上使用它。 - JRP

Answer 6

首先进行线性回归，例如通过普通最小二乘法。这将给出圆柱的轴。

现在计算每个点到该3D轴垂直的所有距离。最大值是圆柱的半径。

如果在计算垂直距离的过程中，您还获得了沿着轴线内部对齐的距离（在某个远离的原点处），那么最小和最大对齐距离就是顶部和底部面的距离。

- Ripi2

---

你有证据证明这是最小值吗？ - Thomas Cohn

不会是最小值：OLS 最小化残差平方和，而你需要最小化最大误差才能解决这个问题。 - c2huc2hu

以下是我的考虑。假设我们已经找到了所需的圆柱体。将3D点集分为两个子集：子集1 =靠近圆柱体表面的薄层中的点，子集2 =在圆柱体内部且比表面更靠近轴线的点。如果子集1更加密集，则您的方法将更精确地逼近轴线位置。但是，如果子集2更加密集，则OLS找到的轴线位置将在圆柱体内任意随机位置。 - JRP

举个反例，假设我有10,000个点在(0,0,0)，10,000个点在(1,0,0)，还有1个点在(0,5,0)。线性回归会将圆柱的轴放得非常靠近x轴，但y轴会产生半径较小的圆柱。 - idontseethepoint