对于光学扫描设备(如激光雷达)采集到的非规则点云数据，最重要的需求之一就是进行表面重建(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->loc);
int number=1;
vtkIdType neighborIndex;
//把所有点坐标加起来除以数量
for(j=0;j<p->neighbors->GetNumberOfIds();j++)
{
neighborIndex = p->neighbors->GetId(j);
pointi = input->GetPoint(neighborIndex);
vtkAddBToA(p->o,pointi);
number++;
}
vtkDivideBy(p->o,number);
// then compute the covariance matrix
//所有点坐标减中心点坐标 求协方差矩阵 加起来
vtkSRMakeZero(covar,0,2,0,2);
for(k=0;k<3;k++)
v3d[k] = p->loc[k] - p->o[k];
vtkSRAddOuterProduct(covar,v3d);
for(j=0;j<p->neighbors->GetNumberOfIds();j++)
{
neighborIndex = p->neighbors->GetId(j);
pointi = input->GetPoint(neighborIndex);
for(k=0;k<3;k++)
{
v3d[k] = pointi[k] - p->o[k];
}
vtkSRAddOuterProduct(covar,v3d);
}
//除以数量
vtkSRMultiply(covar,1.0/number,0,2,0,2);
// then extract the third eigenvector
vtkMath::Jacobi(covar,eigenvalues,eigenvectors);
// third eigenvector (column 2, ordered by eigenvalue magnitude) is plane normal
for(k=0;k<3;k++)
{
p->n[k] = eigenvectors[k][2];
}
}
vtkSRFreeMatrix(covar,0,2,0,2);
vtkSRFreeVector(v3d,0,2);
vtkSRFreeVector(eigenvalues,0,2);
vtkSRFreeMatrix(eigenvectors,0,2,0,2);
}
//--------------------------------------------------------------------------
// 3a. Compute a cost between every pair of neighbors for the MST
// 初步确定相邻点之间的法线关系 作为最小生成树的结点距离
// --------------------------------------------------------------------------
// cost = 1 - |normal1.normal2|
// ie. cost is 0 if planes are parallel, 1 if orthogonal (least parallel)
for(i=0;i<COUNT;i++)
{
SurfacePoint *p = &surfacePoints[i];
p->costs = new double[p->neighbors->GetNumberOfIds()];
// compute cost between all its neighbors
// (bit inefficient to do this for every point, as cost is symmetric)
for(j=0;j<p->neighbors->GetNumberOfIds();j++)
{
p->costs[j] = 1.0 -
fabs(vtkMath::Dot(p->n,surfacePoints[p->neighbors->GetId(j)].n));
}
}
// --------------------------------------------------------------------------
// 3b. Ensure consistency in plane direction between neighbors
// --------------------------------------------------------------------------
// 确定法线朝向的问题与图的最大割问题等价，是NP难问题，因而不大可能求得精确解。使用最小生成树求法线的近似方向。从一点开始，贪心地找方向最接近的结点设为相同方向。
// method: guess first one, then walk through tree along most-parallel
// neighbors MST, flipping the new normal if inconsistent
// to walk minimal spanning tree, keep record of vertices visited and list
// of those near to any visited point but not themselves visited. Start
// with just one vertex as visited. Pick the vertex in the neighbors list
// that has the lowest cost connection with a visited vertex. Record this
// vertex as visited, add any new neighbors to the neighbors list.
vtkIdList *nearby = vtkIdList::New(); // list of nearby, unvisited points
// start with some vertex
int first=0; // index of starting vertex
surfacePoints[first].isVisited = 1;
// add all the neighbors of the starting vertex into nearby
for(j=0;j<surfacePoints[first].neighbors->GetNumberOfIds();j++)
{
nearby->InsertNextId(surfacePoints[first].neighbors->GetId(j));
}
double cost,lowestCost;
int cheapestNearby = 0, connectedVisited = 0;
//邻域待访问结点 repeat until nearby is empty
while(nearby->GetNumberOfIds()>0)
{
// for each nearby point:
vtkIdType iNearby,iNeighbor;
lowestCost = VTK_DOUBLE_MAX;
//对于每个未访问结点
for(i=0;i<nearby->GetNumberOfIds();i++)
{
iNearby = nearby->GetId(i);
// 遍历其邻域的每个已访问结点，计算方向的相似程度 找最小的
for(j=0;j<surfacePoints[iNearby].neighbors->GetNumberOfIds();j++)
{
iNeighbor = surfacePoints[iNearby].neighbors->GetId(j);
if(surfacePoints[iNeighbor].isVisited)
{
cost = surfacePoints[iNearby].costs[j];
// pick lowest cost for this nearby point
if(cost<lowestCost)
{
lowestCost = cost;
cheapestNearby = iNearby;
connectedVisited = iNeighbor;
// 如果基本平行的话不用继续了 直接跳出
if(lowestCost<0.1)
{
i = nearby->GetNumberOfIds();
break;
}
}
}
}
}
//已访问结点=未访问结点 不可能的
if(connectedVisited == cheapestNearby)
{
vtkErrorMacro (<< "Internal error in vtkSurfaceReconstructionFilter");
return 0;
}
// 如果方向相反 就反向
if(vtkMath::Dot(surfacePoints[cheapestNearby].n,
surfacePoints[connectedVisited].n)<0.0F)
{
// flip this normal
vtkMultiplyBy(surfacePoints[cheapestNearby].n,-1);
}
// add this nearby point to visited
surfacePoints[cheapestNearby].isVisited = 1;
// remove from nearby
nearby->DeleteId(cheapestNearby);
// add all new nearby points to nearby 继续在该点的邻域搜索
for(j=0;j<surfacePoints[cheapestNearby].neighbors->GetNumberOfIds();j++)
{
iNeighbor = surfacePoints[cheapestNearby].neighbors->GetId(j);
if(surfacePoints[iNeighbor].isVisited == 0)
{
nearby->InsertUniqueId(iNeighbor);
}
}
}
nearby->Delete();
// --------------------------------------------------------------------------
// 4. Compute signed distance to surface for every point on a 3D grid
// --------------------------------------------------------------------------
{
// 计算输出网格的长宽高、密度
double bounds[6];
for(i=0;i<3;i++)
{
bounds[i*2]=input->GetBounds()[i*2];
bounds[i*2+1]=input->GetBounds()[i*2+1];
}
// estimate the spacing if required
if(this->SampleSpacing<=0.0)
{
// spacing guessed as cube root of (volume divided by number of points)
//体积除以点数 开立方根= 平均每个点的边长
this->SampleSpacing = pow(static_cast<double>(bounds[1]-bounds[0])*
(bounds[3]-bounds[2])*(bounds[5]-bounds[4]) /
static_cast<double>(COUNT),
static_cast<double>(1.0/3.0));
}
// allow a border around the volume to allow sampling around the extremes
for(i=0;i<3;i++)
{
bounds[i*2]-=this->SampleSpacing*2;
bounds[i*2+1]+=this->SampleSpacing*2;
}
double topleft[3] = {bounds[0],bounds[2],bounds[4]};
double bottomright[3] = {bounds[1],bounds[3],bounds[5]};
int dim[3];
for(i=0;i<3;i++)
{
dim[i] = static_cast<int>((bottomright[i]-topleft[i])/this->SampleSpacing);
}
// initialise the output volume
//输出的是全部内容
outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1);
//输出的范围
output->SetExtent(0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1);
//原点坐标
output->SetOrigin(bounds[0], bounds[2], bounds[4]); // these bounds take into account the extra border space introduced above
//沿三个方向的坐标间距
output->SetSpacing(this->SampleSpacing, this->SampleSpacing, this->SampleSpacing);
//根据刚刚设置的参数分配存储空间 由outInfo确定数据类型和通道数
output->AllocateScalars(outInfo);
outInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),
0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1);
vtkFloatArray *newScalars =
vtkFloatArray::SafeDownCast(output->GetPointData()->GetScalars());
outInfo->Set(vtkDataObject::SPACING(),
this->SampleSpacing, this->SampleSpacing, this->SampleSpacing);
outInfo->Set(vtkDataObject::ORIGIN(),topleft,3);
// initialise the point locator (have to use point insertion because we need to set our own bounds, slightly larger than the dataset to allow for sampling around the edge)
//邻域搜索器
vtkPointLocator *locator = vtkPointLocator::New();
vtkPoints *newPts = vtkPoints::New();
//指定数据存储的位置和位置范围盒
locator->InitPointInsertion(newPts,bounds,static_cast<int>(COUNT));
for(i=0;i<COUNT;i++)
{
locator->InsertPoint(i,surfacePoints[i].loc);
}
// go through the array probing the values
int x,y,z;
int iClosestPoint;
int zOffset,yOffset,offset;
double probeValue;
double point[3],temp[3];
for(z=0;z<dim[2];z++)
{
zOffset = z*dim[1]*dim[0];
point[2] = topleft[2] + z*this->SampleSpacing;
for(y=0;y<dim[1];y++)
{
yOffset = y*dim[0] + zOffset;
point[1] = topleft[1] + y*this->SampleSpacing;
for(x=0;x<dim[0];x++)
{
point[0] = topleft[0] + x*this->SampleSpacing;
// find the distance from the probe to the plane of the nearest point
iClosestPoint = locator->FindClosestInsertedPoint(point);
vtkCopyBToA(temp,point);
//方格顶点坐标减去数据点坐标
vtkSubtractBFromA(temp,surfacePoints[iClosestPoint].loc);
//在法线上的投影长度 决定signed distance
probeValue = vtkMath::Dot(temp,surfacePoints[iClosestPoint].n);
offset = x + yOffset; //点索引
newScalars->SetValue(offset,probeValue);
}
}
}
locator->Delete();
newPts->Delete();
}
delete [] surfacePoints;
return 1;
}
Hugues Hoppe在其1994年的博士论文中对于曲面的优化和保角简化还有一整套详细讨论。
vtkExtractSurface使用了类似的方法，并可以利用Range的可见性信息以确定法线方向。
