A quick and rather easy to implement algorithm for getting a decent hull of triangles around the points (outside and inside!) is based on the "Marching Cubes" algorithm:
- define a grid (minimal x,y,z values + grid resolutions in x,y,z direction + number of points in x,y,z direction)
- initialize a value for each grid point to zero
- "rasterize" the given points (round the coordinates to the nearest grid point); for the trivial version just set the corresponding grid point's value to 1
- Now you can polygonize each basic cube of the grid with the "Marching Cubes" algorithm with an
isolevel
of0.9
.
These code snippets might help you. For initializing:
union XYZ {
struct {
double x,y,z;
};
double coord[3];
};
typedef std::vector<XYZ> TPoints;
TPoints points;
// get points from file or another data structure
// initialize these values for your needs (the bounding box of the points will help)
XYZ coordsMin, coordsMax;
XYZ voxelSize; // grid resolution in each coordinate direction
int gridSize; // number of grid points, here the same for each coordinate direction
int gridSize2 = gridSize*gridSize;
char* gridVals = new char[gridSize2*gridSize];
for (int i=0; i<gridSize; ++i)
for (int j=0; j<gridSize; ++j)
for (int k=0; k<gridSize; ++k)
gridVals[i*gridSize2+j*gridSize+k] = 0;
for (size_t i=0; i<points,size(); ++i)
{
XYZ const& p = points[i];
int gridCoords[3];
for (int c=0; c<3; ++c)
gridCoords[c] = (int)((p.coord[c]-coordsMin.coord[c])/voxelSize.coord[c]);
int gridIdx = gridCoords[0]*gridSize2 + gridCoords[1]*gridSize + gridCoords[2];
gridVals[gridIdx] = 1;
}
And then for computing the triangles based on the function Polygonize
from the cited implementation:
const double isolevel = 0.9;
TRIANGLE triangles[5]; // maximally 5 triangles for each voxel
int cellCnt = 0;
for (int i=0; exportOk && i<gridSize-1; ++i)
for (int j=0; exportOk && j<gridSize-1; ++j)
for (int k=0; exportOk && k<gridSize-1; ++k)
{
GRIDCELL cell;
for (int cornerIdx=0; cornerIdx<8; ++cornerIdx)
{
XYZ& corner = cell.p[cornerIdx];
// the function Polygonize expects this strange order of corner indexing ...
// (see http://paulbourke.net/geometry/polygonise/)
int xoff = ((cornerIdx+1)/2)%2;
int yoff = (cornerIdx/2)%2;
int zoff = cornerIdx/4;
corner.x = (i+xoff)*voxelSize.x + coordsMin.x;
corner.y = (j+yoff)*voxelSize.y + coordsMin.y;
corner.z = (k+zoff)*voxelSize.z + coordsMin.z;
cell.val[cornerIdx] = gridVals[(i+xoff)*gridSize2+(j+yoff)*gridSize+k+zoff];
}
int triangCnt = Polygonise(cell, isolevel, triangles);
triangCntTotal += triangCnt;
triangCntTotal += triangCnt;
for (int t=0; t<triangCnt; ++t)
{
TTriangle const& triangle = triangles[t].p;
ExportTriangle(triangle);
}
}
This did the trick for me in an industrial application.