vtk中的点云曲面重建

对于光学扫描设备(如激光雷达)采集到的非规则点云数据，最重要的需求之一就是进行表面重建(Surface Reconstruction)：对成片密集分布的点云以三角片拟合，形成连续、精确、良态的曲面表示。
目前主流的算法可分为剖分类、组合类和拟合类。剖分类比如Voronoi图、Delaunay三角剖分，原始数据点即为顶点，数据无损失，数据冗余多，生成的曲面不光顺，容易受噪声影响。组合类比如 Power Crust、DBRG在剖分的基础上使用简单几何形状组合，使算法更稳定，并提供了拓扑信息。拟合类的算法以Hugues Hoppe提出的Signed Distance Function和Poisson方法为代表。所得曲面并不完全与数据点重合，而是进行了一定的合并和简化，并对噪声有一定鲁棒性。
经典的Signed Distance Function重建算法主要流程如下：

1. 对每个数据点，搜索其邻域数据点，使用特征向量方法计算法线和法平面
2. 由于法线的方向可有正负两个不确定，故对全局使用最小生成树近似计算法线朝向
3. 以立方体素(voxel)为单位，沿着曲面滑动延展，计算每个格顶点到法平面的距离
4. 使用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->