对于光学扫描设备(如激光雷达)采集到的非规则点云数据,最重要的需求之一就是进行表面重建(Surface Reconstruction):对成片密集分布的点云以三角片拟合,形成连续、精确、良态的曲面表示。
目前主流的算法可分为剖分类、组合类和拟合类。剖分类比如Voronoi图、Delaunay三角剖分,原始数据点即为顶点,数据无损失,数据冗余多,生成的曲面不光顺,容易受噪声影响。组合类比如 Power Crust、DBRG在剖分的基础上使用简单几何形状组合,使算法更稳定,并提供了拓扑信息。拟合类的算法以Hugues Hoppe提出的Signed Distance Function和Poisson方法为代表。所得曲面并不完全与数据点重合,而是进行了一定的合并和简化,并对噪声有一定鲁棒性。
经典的Signed Distance Function重建算法主要流程如下:
- 对每个数据点,搜索其邻域数据点,使用特征向量方法计算法线和法平面
- 由于法线的方向可有正负两个不确定,故对全局使用最小生成树近似计算法线朝向
- 以立方体素(voxel)为单位,沿着曲面滑动延展,计算每个格顶点到法平面的距离
- 使用Marching Cube算法从立方体素中(插值地)提取三角面
VTK库提供的vtkSurfaceReconstructionFilter类实现了该算法。下面对其源码进行注释
int vtkSurfaceReconstructionFilter::RequestData(
vtkInformation* vtkNotUsed( request ),
vtkInformationVector** inputVector,
vtkInformationVector* outputVector)
{
// 从pipeline中获取输入接口
vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
vtkDataSet *input = vtkDataSet::SafeDownCast(
inInfo->Get(vtkDataObject::DATA_OBJECT()));
// 从pipeline中获取输出接口
vtkInformation *outInfo = outputVector->GetInformationObject(0);
vtkImageData *output = vtkImageData::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
const vtkIdType COUNT = input->GetNumberOfPoints();
SurfacePoint *surfacePoints;
surfacePoints = new SurfacePoint[COUNT];
vtkIdType i, j;
int k;
// --------------------------------------------------------------------------
// 1. Build local connectivity graph
// -------------------------------------------------------------------------
{
//八叉树 用于邻域搜索
vtkPointLocator *locator = vtkPointLocator::New();
locator->SetDataSet(input);
vtkIdList *locals = vtkIdList::New();
// if a pair is close, add each one as a neighbor of the other
for(i=0;i<COUNT;i++)
{
//遍历所有点
SurfacePoint *p = &surfacePoints[i];
vtkCopyBToA(p->loc,input->GetPoint(i));
//查找当前点的邻域
locator->FindClosestNPoints(this->NeighborhoodSize,p->loc,locals);
int iNeighbor;
for(j=0;j<locals->GetNumberOfIds();j++)
{
iNeighbor = locals->GetId(j);
if(iNeighbor!=i)
{
//邻域有点 双向记录
p->neighbors->InsertNextId(iNeighbor);
surfacePoints[iNeighbor].neighbors->InsertNextId(i);
}
}
}
locator->Delete();
locals->Delete();
}
// --------------------------------------------------------------------------
// 2. Estimate a plane at each point using local points
// --------------------------------------------------------------------------
{
double *pointi;
double **covar,*v3d,*eigenvalues,**eigenvectors;
covar = vtkSRMatrix(0,2,0,2);
v3d = vtkSRVector(0,2);
eigenvalues = vtkSRVector(0,2);
eigenvectors = vtkSRMatrix(0,2,0,2);
for(i=0;i<COUNT;i++)
{
SurfacePoint *p = &surfacePoints[i];
// first find the centroid of the neighbors
vtkCopyBToA(p->o,p->