Reputation: 113
Expand 3d convex polygon
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.
Upvotes: 3
Views: 434
Answers (1)
Reputation: 51863
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,i2of 3 adjacent non parallel triangles to each point. I also addedmidpoint of each triangle anddnormal 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,midfor each face that takesO(n)assumingntriangles andmpoints. And the adjacency indexes for each point are computed inO(m)so the whole thing isO(n+m)The adjacency is computed easily first cleari0,i1,i2of 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
midpoint bynormal*offset_steprecomputingdfor all faces. After that you loop through all points and compute intersection of 3 planes you got index to. So this is alsoO(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)
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:
Upvotes: 1
Related Questions
- Minkowski Sum of 3D(+) convex polygons
- Create a 3D polyhedron from 2D convex polygon in bullet physics
- Expand polygons with boost::geometry?
- Expand fill of convex polygon
- Expanding Convex Hull to Reduce Sides
- Simplified concave-hulls
- Polygon math
- Representing a convex hull in higher dimensions (3+)
- Convex 3D polygonal object
- Distance between two convex polygons in 3D