Ответ 1

Сначала найдите ориентированную ограничительную шкатулу (OBB) для облака точек. Для этого есть несколько алгоритмов. Вероятно, этот вариант оптимален:

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

Теперь неоптимальный ориентированный цилиндр, охватывающий OBB, можно легко найти, вращая OBB вокруг его самой длинной оси. Аналогично, цилиндр, заключенный OBB, может быть найден как имеющий ту же ось, что и другой, но радиус равен половине кратчайшей стороны OBB, нормальной к оси.

Моя цель заключается в том, что оптимальный радиус цилиндра находится между этими двумя цилиндрами.

Лучший цилиндр можно легко найти, если вы вычислили минимальное расстояние всех точек до внешнего цилиндра и отрегулировали его радиус, чтобы сделать это расстояние равным нулю.

Этот подход, вероятно, работает, но не является оптимально вычислительным, поскольку вам нужно вычислить расстояние от всех точек до цилиндра. Возможно, внутренний цилиндр может использоваться для отбраковки всех точек, находящихся внутри него. Я не слишком много думал об этой идее.

ОБНОВИТЬ:

Похоже, что вопрос не ясен о том, что является "наименьшим" и на самом деле требует вещей за "наименьшим", и это не очень хорошо. Предполагается, что "самый маленький" цилиндр, охватывающий облако точек, минимизирует пустое пространство внутри цилиндра (по крайней мере, я понимаю это как наименьшее). Но ОП также накладывает ограничение на то, что наименьший цилиндр должен соответствовать форме входных данных. Это означает, что если входные данные являются выборкой на половину цилиндра (разрезанной самой одинокой стороной), ответом должен быть цилиндр, который лучше всего соответствует форме этой половины. Независимо от того, имеет ли этот цилиндр более пустое пространство, чем другой цилиндр, который охватывает данные.

Эти два требования являются противоречивыми. Так как наименьший цилиндр может не соответствовать изогнутой форме данных и лучше всего подходит цилиндр, изогнутая форма данных не может быть наименьшим цилиндром.

Мой ответ (и другие ответы), основанный на OBB, отвечает на вопрос относительно "самого маленького" цилиндра, охватывающего данные, минимизирующие пустое пространство внутри цилиндра.

В другом случае при установке цилиндра в форму данных можно также ответить, используя подход оптимизации. Но ответа нет. "Лучший" цилиндр зависит от потребностей приложения и должен вычисляться по крайней мере с двумя разными стратегиями в зависимости от данных.

Ответ 2

1. вычислить OBB

   поэтому либо используйте СПС, либо это

   • Как вычислить OBB нескольких кривых?

   для получения 3D OBB. Код в ссылке должен быть перенесен в 3D, но принцип тот же.

2. первоначальное предположение

   поэтому OBB предоставит вам ориентированную ограничительную рамку. Его самая большая сторона будет параллельна оси вращения вашего цилиндра. Итак, давайте начнем с цилиндровым outscribing это О. Таким образом, центральная ось будет центром OBB и параллельна ее самой большой стороне. (если у вас нет самой большой стороны, вам нужно проверить все 3 комбинации). Диаметр будет больше по сравнению с оставшимися сторонами.

3. подходящий цилиндр

   Теперь попробуйте "все" комбинации смещения и радиуса (может быть также высоты), охватывающие все ваши точки вблизи первоначальной догадки и запомните лучший (в соответствии с вашими желаемыми спецификациями). Вы можете использовать любой метод оптимизации для этого, но мой любимый:

   • Как работает поиск приближений

Действительность результата зависит от процесса подгонки. Но не становитесь слишком дикими с вложенными фитингами, так как сложность тоже дикая.

[Edit1] 3D OBB в C++

Мне было любопытно, и я получил некоторое время сегодня, поэтому я закодировал 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[] на каждую ось или в виде 8 угловых точек OBB p

Материал работает просто, используя "все" сферические положения с некоторым шагом для оси w а затем попробуйте все полярные положения u,v перпендикулярные каждому, и w помнящие минимальный объем OBB. Затем рекурсивно искать только положения, найденные наилучшим решением, с меньшим шагом для повышения точности.

Я думаю, что это должно обеспечить прекрасную начальную точку. Если вы используете минимальный круг вместо u,v вращение (цикл for (ec=1,c=c0;ec;c+=dc)), то вы можете получить свой цилиндр непосредственно из этого поиска.

Код еще не оптимизирован (некоторые части, такие как проверка оси w) можно переместить на нижний уровень вложенного цикла. Но я хотел сохранить это простым и понятным, насколько мог, вместо этого.

[Edit2] 3D OBC в C++

Мне удалось изменить свой 3D OBB, заменив U,V поиск минимальным окружением (надеюсь, что я правильно его реализовал, но это похоже на это...), которые находят минимальный окружающий 2D-круг всех точек, проектируемых на UV плоскости, что делает его ориентированный ограничивающий цилиндр, параллельный W Я использовал первый подход из PDF из вашей ссылки (используя биссектрису). Вот результат:

В синем - 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);
    }

Where List<double> pnt - мой шаблон динамического массива double pnt[]. Что здесь не важно.

Помните, что если вы выбрали слишком большой начальный шаг (da,db) для поиска направления W вы можете пропустить правильное решение, заманив себя внутри локального минимума.

Ответ 3

Найти точное решение представляется очень сложной задачей. Сделав гипотезы о направлении оси и вращая облако (из которого вы держите только вершины выпуклой оболочки) и проецируя точки на XY, вы действительно можете обратиться к минимальной проблеме окружного круга.

Но вам нужно сделать несколько гипотез о возможных направлениях. Для точного решения вы должны попробовать все направления, чтобы окружность окружения определялась разными парами или тройками вершин. Это определяет сложные области в пространстве углов поворота, и для каждой области есть точка, которая достигает оптимального. Это связано с проблемой нелинейной минимизации с нелинейными ограничениями и кажется едва выполнимой.

На этом этапе все, что я могу порекомендовать, является приблизительным решением, так что вы пытаетесь фиксированное количество предопределенных направлений и решаете соответствующий круг. Так как решение в любом случае приблизительное, может применяться приблизительная окружность.

Ответ 4

Концептуальный ответ

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 используемый 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

Наименьший объем цилиндра, который может содержать эти точки,

Заметки

• Максимальное расстояние между точками было найдено при использовании функции pdist. Затем он был отформатирован с squareform чтобы выбросить его в Pandas DataFrame поэтому индексы обеих точек можно легко найти с помощью idxmax. Вероятно, это лучший способ сделать это без Панд.

• Если часть np.cross оставила вас почесывать голову, это просто найти минимальное расстояние между точкой и линией. Я могу следить за деталями, если вам интересно, но если вы рисуете перекрестное произведение двух строк, вы получаете параллелограмм, где две непараллельные стороны задаются линиями. Этот параллелограмм имеет ту же площадь, что и прямоугольник с длиной, равной одной из линий, и шириной, равной расстоянию от точки до линии.

Ответ 5

Алгоритм грубой силы

Начнем с самого простого случая: 3 балла.

Самый маленький цилиндр имеет три точки на своей поверхности. Наименьший радиус означает, что две точки находятся в диаметре поперечного сечения, даже если эта секция не перпендикулярна оси цилиндра. То же самое касается третьего пункта.
Таким образом, ось проходит через центр диаметра, который является срединной точкой сегмента, определяемой двумя точками.
Таким образом, ось определяется двумя средними точками.

У нас есть три средних точки и три возможные оси:

Наилучшее решение (минимальный радиус) выбирается путем нахождения минимума ri среди решений.

Заметим, что каждая ось ai параллельна некоторому сегменту Pi,Pj.

GENERIC, n POINTS CASE

Скажем, мы нашли решение. По крайней мере три точки этого решения лежат на поверхности, что аналогично случаю с тремя точками. Если мы найдем триплет, то мы можем применить метод средних точек.
Таким образом, ось решения представляет собой некоторую линию через две средние точки. В этом алгоритме грубой силы мы проверяем их все. Для каждой оси вычислите расстояние всех точек, перпендикулярных этой оси, и получим наибольшее значение di, поэтому мы получим радиус окружения для каждой оси.

Решение - это ось для минимума di. Радиус - это минимум di.

Каков порядок этого алгоритма?
Для n точек C(n,2)= n(n-1)/2 средних точек. И ось n(n-1)(n(n-1)-2)/8. Для каждой оси мы проверяем n радиус. Таким образом, мы получаем O (n ⁵)

улучшения

• Постройте сначала выпуклый корпус и отбросьте все внутренние точки
• У меня такое ощущение, что ось для решения определяется самым длинным сегментом Mi,Mj.

РЕДАКТИРОВАТЬ

Как показал @JRP, этот метод не находит оптимального решения для призмы. Попытка с центрами треугольников (вместо середины сегментов) тоже не работает, подумайте о облаке точек с 8 точками в вершинах куба.

Это может быть хорошим аппроксимацией, но не оптимальным решением.

Ответ 6

Бесплатный CYL может сделать это: http://petitjeanmichel.free.fr/itoweb.petitjean.freeware.html; метод описан в моей статье: Об алгебраических решениях наименьших задач вмещающих цилиндров. Appl. Alg. Eng. Comm. Comp., 2012, 23 [3-4], 151-164. doi: 10.1007/s00200-012-0171-y (препринт в arXiv: 1008.5259). Время вычислений может быть большим. Рекомендуется удалить точки внутри выпуклой оболочки, прежде чем вычислять наименьший вмещающий цилиндр.

Ответ 7

Сначала сделайте линейную регрессию, например, по обычным наименьшим квадратам. Это даст вам ось цилиндра.

Теперь вычислите все расстояния от каждой точки, перпендикулярной этой трехмерной оси. Максимальное значение - радиус цилиндра.

Если в процессе вычисления перпендикулярных расстояний вы также выровняете расстояния внутри оси (поместите источник где-то далеко), то минимальные и максимальные выровненные расстояния будут равны для верхней и нижней поверхностей.