Clever Geek Handbook

Expand 3d Convex Polygon

I have a convex polyhedron made up of points and triangles (the normal of the triangle is outside the polygon):

enter image description here

Input structure:

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

I would like to expand the polygon by moving the triangles along their normal some distance. How can I calculate new 3D coordinates of the moved points with good performance?

=> I have an implementation in O (n ^ 3): I move all planes (triangles) along their normal and then I find all points by checking all plane intersections (three planes are not parallel, give an intersection point) This solution works but gives a very poor result.

=> I am also trying to consistently move the point to all normals of its attached triangles, but not giving the correct result.

+3


source to share


1 answer


I was curious, so I tried to describe it in C ++ here:

  • <strong> structures

    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; }
    };

    So, I added indices i0,i1,i2

    of 3 adjacent non-parallel triangles to each point. I've also added the point of mid

    each triangle and the d

    normal offset to speed up the calculations, but both can be calculated when needed, so you don't really need to add them to the grid itself.

  • pre-computation

    so you need to pre-compute nor,d,mid

    for each face receiving O(n)

    , assuming triangles n

    and m

    . And the adjacency indices for each point are calculated in O(m)

    , so that's all O(n+m)

    . Adjacency is easy to calculate at first clear i0,i1,i2

    all points. Then draw on all faces and for each add its index to each of its points, if there are less than 3 normals and the normal is not parallel to it.

  • bias

    now the displacement is performed only by shifting the point mid

    by normal*offset_step

    recalculating d

    for all faces. After that, you iterate over all the points and compute the intersection of the three planes to which you got the index. So this is also O(n+m)

    .

    I was too lazy to get the intersection equation, so I used a 3x3 inverse matrix instead. Since my matrices are 4x4, the last row and column are not used. Beware that my matrices are OpenGL so they are transposed ... which is why the normals are loaded so strangely into it.

Here's my C ++ source:

//---------------------------------------------------------------------------
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);
            }
        }
    };
//---------------------------------------------------------------------------

Here's a preview (left is source and then offset applied multiple times)

preview



Works with negative steps too. Using this looks like this:

// 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);

I am also using my dynamic list template, so:


List<double> xxx;

matches double xxx[];


xxx.add(5);

add 5

to the end of the list
xxx[7]

an accessor array element (safe)
xxx.dat[7]

accessor array element (unsafe but fast direct access)
xxx.num

- the actual used array size
xxx.reset()

clears the array and sets xxx.num=0


xxx.allocate(100)

preset space for 100

elements

And the math source of vector and matrix can be found here:

+1


source






More articles:


All Articles