Expand 3d convex polygon

Question

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

The input structure is:

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 triangles along their normal of a certain distance. How can I compute the new 3D coordinates of 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 the points by testing all planes intersections (3 planes not parallel give an intersection point). This solution works but give really bad performance result.

=> I have also try to move the point sequentially to all normals of his attached triangles but it doesn't give the correct result.

Moving the planes results in a new computation of halfspace intersection, which can be computed more quickly than $n^3$. So maybe just implement halfspace intersection. — Joseph O'Rourke, Apr 28 '17 at 21:29
For each Point I would compute a LUT to which triangles it belongs. This way while resizing you can access all the neighbor planes in `O(1)` The construction of such LUT will take `O(n)` where `n` is number of triangles. Also my guts telling me it should be enough to store 3 non coplanar triangles per point (no need to store all neighbors).... — Spektre, Apr 29 '17 at 09:15

score 1 · Accepted Answer · edited May 23 '17 at 12:26

I was curious so I tried to code this in C++ here some insights:

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 indexes i0,i1,i2 of 3 adjacent non parallel triangles to each point. I also added mid point of each triangle and d normal offset to speed up computations but both can be computed when needed so you do not really need them to add to the mesh itself.

pre-computation

so you need to pre-compute nor,d,mid for each face that takes O(n) assuming n triangles and m points. And the adjacency indexes for each point are computed in O(m) so the whole thing is O(n+m) The adjacency is computed easily first clear i0,i1,i2 of all points. Then loop through all faces and for each add its index to each of its points if there are less than 3 normals and no normal is parallel to it.
offset

the offset is now done just by offsetting mid point by normal*offset_step recomputing d for all faces. After that you loop through all points and compute intersection of 3 planes you got index to. So this is also O(n+m).

I was too lazy to derive intersection equation so I used 3x3 inverse matrix instead. As my matrices are 4x4 the last row and column is unused. Beware my matrices are OpenGL like so they are transposed... that is why the normals are loaded so weirdly into it.

Here 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 preview (left is source and then few times applied offset)

preview

Works also with negative steps. Ussage of 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 also use mine dynamic list template so:

List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items

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

Understanding 4x4 homogenous transform matrices

Thanks for great answer. I was wondering if it's not possible to avoid matrix multiplication. Indeed, in 2D, I'm able to compute the new point with only: normal of neighbors(firstNormal & secondNormal) and the origin points (point[i]) with this simple formula: double moveLength = offset / (1.0 + firstNormal.dotProduct(secondNormal)); Point2 newPoint(points[i].X + (firstNormal.X + secondNormal.X)*moveLength, points[i].Y + (firstNormal.Y + secondNormal.Y)*moveLength)); — Gérard Legrand, Apr 30 '17 at 11:31
@GérardLegrand re-enabled preview image (I have no clue how the image got converted to code) also I would not bother with the matrix multiplication `3x3 matrix * 3D vector` is fast... more worried I would be with the `3x3` inverse matrix I am using 4x4 with subdeterminants so it is far from optimal but I am too lazy to rewrite that stuff ... — Spektre, Apr 30 '17 at 13:02

Expand 3d convex polygon

1 Answers1