math - 展开 3d 凸多边形

Question

我有一个由点和三角形组成的凸多面体(三角形的法线在多边形之外):

输入结构为:

class Vector {
  float X, Y, Z;
};
class Triangle {
  int pointIndexes[3];
  Vector normal;
};
class Polyhedron{
  vector<Vector> points;
  vector<Triangle> triangles;
};

我想通过沿着三角形的法线移动一定距离来扩展多边形。如何以良好的性能计算移动点的新 3D 坐标？

=> 我有一个 O(n^3) 实现:我沿着法线移动所有平面(三角形)，然后通过测试所有平面相交找到所有点(3 个不平行的平面给出一个相交点) 。该解决方案有效，但性能结果非常糟糕。

=> 我还尝试将点按顺序移动到其附加三角形的所有法线，但它没有给出正确的结果。

Answer 1

我很好奇，所以我尝试用C++对其进行编码，这里有一些见解:

1. 结构

   struct _pnt
       {
       int i0,i1,i2;   // neighbor triangles (non parallel)
       double pos[3];  // position
       _pnt()      {}
       _pnt(_pnt& a)   { *this=a; }
       ~_pnt() {}
       _pnt* operator = (const _pnt *a) { *this=*a; return this; }
       //_pnt* operator = (const _pnt &a) { ...copy... return this; }
       };
   
   struct _fac
       {
       int i0,i1,i2;       // triangle
       double nor[3];      // normal
       double mid[3],d;    // mid point x,y,z and normal d
       _fac()      {}
       _fac(_fac& a)   { *this=a; }
       ~_fac() {}
       _fac* operator = (const _fac *a) { *this=*a; return this; }
       //_fac* operator = (const _fac &a) { ...copy... return this; }
       };
   

   所以我添加了索引 i0,i1,i2每个点有 3 个相邻的不平行三角形。我还添加了mid每个三角形的点和 d正常偏移以加速计算，但两者都可以在需要时计算，因此您实际上并不需要将它们添加到网格本身。

2. 预计算

   所以你需要预先计算nor,d,mid对于每张需要O(n)的脸假设n三角形和m点。每个点的邻接索引在 O(m) 中计算。所以整个事情是O(n+m)首先清楚邻接关系很容易计算 i0,i1,i2所有点。然后循环遍历所有面，如果法线少于 3 个且没有法线与之平行，则为每个面添加其索引到每个点。

3. 偏移量

   现在只需通过偏移 mid 即可完成偏移点 normal*offset_step重新计算d对于所有面孔。之后，您循环遍历所有点并计算您获得索引的 3 个平面的交集。所以这也是O(n+m) .

   我懒得推导交集方程，所以我用 3x3 逆矩阵代替。由于我的矩阵是 4x4，所以最后一行和最后一列未使用。请注意，我的矩阵与 OpenGL 类似，因此它们会转置...这就是法线加载到其中的原因。

这是我的 C++ 源代码:

//---------------------------------------------------------------------------
struct _pnt
    {
    int i0,i1,i2;   // neighbor triangles (non parallel)
    double pos[3];  // position
    _pnt()      {}
    _pnt(_pnt& a)   { *this=a; }
    ~_pnt() {}
    _pnt* operator = (const _pnt *a) { *this=*a; return this; }
    //_pnt* operator = (const _pnt &a) { ...copy... return this; }
    };
//---------------------------------------------------------------------------
struct _fac
    {
    int i0,i1,i2;       // triangle
    double nor[3];      // normal
    double mid[3],d;    // mid point x,y,z and normal d
    _fac()      {}
    _fac(_fac& a)   { *this=a; }
    ~_fac() {}
    _fac* operator = (const _fac *a) { *this=*a; return this; }
    //_fac* operator = (const _fac &a) { ...copy... return this; }
    };
//---------------------------------------------------------------------------
class mesh
    {
public:
    List<_pnt> pnt;
    List<_fac> fac;

    mesh()      {}
    mesh(mesh& a)   { *this=a; }
    ~mesh() {}
    mesh* operator = (const mesh *a) { *this=*a; return this; }
    //mesh* operator = (const mesh &a) { ...copy... return this; }

    void icosahedron(double r)  // init mesh with icosahedron of radius r
        {
        // points
        double a=r*0.525731112119133606;
        double b=r*0.850650808352039932;
        _pnt p; p.i0=-1; p.i1=-1; p.i2=-1; pnt.num=0;
        vector_ld(p.pos,-a,0.0, b); pnt.add(p);
        vector_ld(p.pos, a,0.0, b); pnt.add(p);
        vector_ld(p.pos,-a,0.0,-b); pnt.add(p);
        vector_ld(p.pos, a,0.0,-b); pnt.add(p);
        vector_ld(p.pos,0.0, b, a); pnt.add(p);
        vector_ld(p.pos,0.0, b,-a); pnt.add(p);
        vector_ld(p.pos,0.0,-b, a); pnt.add(p);
        vector_ld(p.pos,0.0,-b,-a); pnt.add(p);
        vector_ld(p.pos, b, a,0.0); pnt.add(p);
        vector_ld(p.pos,-b, a,0.0); pnt.add(p);
        vector_ld(p.pos, b,-a,0.0); pnt.add(p);
        vector_ld(p.pos,-b,-a,0.0); pnt.add(p);
        // triangles
        _fac f; fac.num=0; vector_ld(f.nor,0.0,0.0,0.0);
        f.i0= 0; f.i1= 4; f.i2= 1; fac.add(f);
        f.i0= 0; f.i1= 9; f.i2= 4; fac.add(f);
        f.i0= 9; f.i1= 5; f.i2= 4; fac.add(f);
        f.i0= 4; f.i1= 5; f.i2= 8; fac.add(f);
        f.i0= 4; f.i1= 8; f.i2= 1; fac.add(f);
        f.i0= 8; f.i1=10; f.i2= 1; fac.add(f);
        f.i0= 8; f.i1= 3; f.i2=10; fac.add(f);
        f.i0= 5; f.i1= 3; f.i2= 8; fac.add(f);
        f.i0= 5; f.i1= 2; f.i2= 3; fac.add(f);
        f.i0= 2; f.i1= 7; f.i2= 3; fac.add(f);
        f.i0= 7; f.i1=10; f.i2= 3; fac.add(f);
        f.i0= 7; f.i1= 6; f.i2=10; fac.add(f);
        f.i0= 7; f.i1=11; f.i2= 6; fac.add(f);
        f.i0=11; f.i1= 0; f.i2= 6; fac.add(f);
        f.i0= 0; f.i1= 1; f.i2= 6; fac.add(f);
        f.i0= 6; f.i1= 1; f.i2=10; fac.add(f);
        f.i0= 9; f.i1= 0; f.i2=11; fac.add(f);
        f.i0= 9; f.i1=11; f.i2= 2; fac.add(f);
        f.i0= 9; f.i1= 2; f.i2= 5; fac.add(f);
        f.i0= 7; f.i1= 2; f.i2=11; fac.add(f);
        // precompute stuff
        compute();
        }
    void compute()  // compute normals and neighbor triangles info
        {
        int i,j,k;
        _fac *f,*ff;
        _pnt *p;
        double a[3],b[3];
        const double nor_dot=0.001;     // min non parallel dor product of normals
        // normals
        for (f=fac.dat,i=0;i<fac.num;i++,f++)
            {
            vector_sub(a,pnt.dat[f->i1].pos,pnt.dat[f->i0].pos);
            vector_sub(b,pnt.dat[f->i2].pos,pnt.dat[f->i0].pos);
            vector_mul(a,b,a);
            vector_one(f->nor,a);
            }
        // mid points
        for (f=fac.dat,i=0;i<fac.num;i++,f++)
            {
            // mid point of triangle as start point
            vector_copy(a,  pnt.dat[f->i0].pos);
            vector_add (a,a,pnt.dat[f->i1].pos);
            vector_add (a,a,pnt.dat[f->i2].pos);
            vector_mul (f->mid,a,0.33333333333);
            f->d=vector_mul(f->mid,f->nor);
            }
        // neighbors
        for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
            {
            p->i0=-1;
            p->i1=-1;
            p->i2=-1;
            }
        for (f=fac.dat,i=0;i<fac.num;i++,f++)
            {
            for (p=pnt.dat+f->i0;p->i2<0;)
                {
                if (p->i0>=0) { ff=fac.dat+p->i0; if (fabs(vector_mul(f->nor,ff->nor))<=nor_dot) break; } else { p->i0=i; break; }
                if (p->i1>=0) { ff=fac.dat+p->i1; if (fabs(vector_mul(f->nor,ff->nor))<=nor_dot) break; } else { p->i1=i; break; }
                                                                                                                 p->i2=i; break;
                }
            for (p=pnt.dat+f->i1;p->i2<0;)
                {
                if (p->i0>=0) { ff=fac.dat+p->i0; if (fabs(vector_mul(f->nor,ff->nor))<=nor_dot) break; } else { p->i0=i; break; }
                if (p->i1>=0) { ff=fac.dat+p->i1; if (fabs(vector_mul(f->nor,ff->nor))<=nor_dot) break; } else { p->i1=i; break; }
                                                                                                                 p->i2=i; break;
                }
            for (p=pnt.dat+f->i2;p->i2<0;)
                {
                if (p->i0>=0) { ff=fac.dat+p->i0; if (fabs(vector_mul(f->nor,ff->nor))<=nor_dot) break; } else { p->i0=i; break; }
                if (p->i1>=0) { ff=fac.dat+p->i1; if (fabs(vector_mul(f->nor,ff->nor))<=nor_dot) break; } else { p->i1=i; break; }
                                                                                                                 p->i2=i; break;
                }
            }
        }
    void draw() // render mesh
        {
        int i;
        _fac *f;
        glBegin(GL_TRIANGLES);
        for (f=fac.dat,i=0;i<fac.num;i++,f++)
            {
            glNormal3dv(f->nor);
            glVertex3dv(pnt.dat[f->i0].pos);
            glVertex3dv(pnt.dat[f->i1].pos);
            glVertex3dv(pnt.dat[f->i2].pos);
            }
        glEnd();
        }
    void draw_normals(double r) // render mesh normals as line of size r
        {
        int i;
        double a[3];
        _fac *f;
        for (f=fac.dat,i=0;i<fac.num;i++,f++)
            {
            // normal endpoints
            vector_mul (a,r,f->nor);
            vector_add (a,a,f->mid);
            glBegin(GL_LINES);
            glVertex3dv(f->mid);
            glVertex3dv(a);
            glEnd();
            // wireframe
            glBegin(GL_LINE_LOOP);
            glVertex3dv(pnt.dat[f->i0].pos);
            glVertex3dv(pnt.dat[f->i1].pos);
            glVertex3dv(pnt.dat[f->i2].pos);
            glEnd();
            }
        }
    void offset(double d)   // offset triangles by d in normal direction
        {
        int i,j;
        _fac *f;
        _pnt *p;
        double a[3],m[16];
        // translate mid points
        for (f=fac.dat,i=0;i<fac.num;i++,f++)
            {
            vector_mul(a,d,f->nor);
            vector_add(f->mid,f->mid,a);
            f->d=vector_mul(f->mid,f->nor);
            }
        // recompute points as intersection of 3 planes
        for (p=pnt.dat,i=0;i<pnt.num;i++,p++)
         if (p->i2>=0)
            {
            matrix_one(m);
            for (f=fac.dat+p->i0,j=0;j<3;j++) m[0+(j<<2)]=f->nor[j]; a[0]=f->d;
            for (f=fac.dat+p->i1,j=0;j<3;j++) m[1+(j<<2)]=f->nor[j]; a[1]=f->d;
            for (f=fac.dat+p->i2,j=0;j<3;j++) m[2+(j<<2)]=f->nor[j]; a[2]=f->d;
            matrix_inv(m,m);
            matrix_mul_vector(p->pos,m,a);
            }
        }
    };
//---------------------------------------------------------------------------

此处预览(左侧是源，然后是几次应用的偏移)

也适用于负面步骤。它的用法如下所示:

// globals
mesh obj;
// init
obj.icosahedron(0.5);
// render   
glFrontFace(GL_CW); // for gluPerspective
glEnable(GL_CULL_FACE);
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
glEnable(GL_COLOR_MATERIAL);
glColor3f(0.5,0.5,0.5);
obj.draw();
glDisable(GL_LIGHTING);
glColor3f(0.0,0.9,0.9);
glLineWidth(2.0);
obj.draw_normals(0.2);
glLineWidth(1.0);
// on mouse wheel
if (WheelDelta>0) obj.offset(+0.05);
if (WheelDelta<0) obj.offset(-0.05);

我也使用我的动态列表模板:

List<double> xxx;与 double xxx[]; 相同
xxx.add(5);添加5到列表末尾
xxx[7]访问数组元素(安全)
xxx.dat[7]访问数组元素(不安全但快速直接访问)
xxx.num是数组实际使用的大小
xxx.reset()清除数组并设置 xxx.num=0
xxx.allocate(100)为 100 预分配空间项目

向量和矩阵数学源可以在这里找到:

• Understanding 4x4 homogenous transform matrices