计算机视觉
图像处理

图形处理(九)点云重建(下)法矢求取、有向距离场等值面提取

一、相关理论

点云重建算法中,点云法矢的求取是非常重要的一步。之前导师让我做点云尖锐特征边重建时,发了一堆异向法矢求取的英文paper,当时我很迷糊,就问了老师,让我做点云重建,为什么发的文献资料都是关于法矢求取的。原来法矢求取,直接关系着重建效果,废话不多说。

点云重建的一般步骤是:

1、点云预处理  。这一步主要是点云滤波、采样、移除离群点等操作,涉及到的经典算法有MLS、双边滤波、WLOP等,这些在我的另外一篇博文中,有比较详细的介绍。

 2、重建网格曲面。这一步主要涉及相关经典算法是点云法矢求取PCA、有向距离场、隐函数、MC算法等相关概念。这篇博文将详解这些算法的实习过程。

二、重建算法流程

重建网格曲面算法流程:

1、通过KD树,求取点云的每个点Pi的k近邻点,对于k值得选取,需要根据点云的密度进行自适应效果会比较好。我比较常用的值是8、16……

这一步可以从网上下载flann库,进行调用,当然如果你熟悉PCL点云开源处理库,那后续的步骤都只是调用一个函数的事了。这一步的代码就不贴了,就是把点云的每个点的k近邻求取出来就完事了。

2、根据k近邻进行计算点云法矢。法矢的初步估算这一步经典的算法是PCA算法,PCA算法是点云法矢计算的常用算法。

a.均匀权的pi点法矢计算公式为:


其中pj为p点的k近邻。然后通过求解协方差矩阵M的最小特征值的特征向量,作为p点的法向量,这就是所谓的主元分析方法PCA。

b.基于高斯权重的主元分析方法:

原来的方法,采用的是均匀权重的方法(1/k),把均匀权改为高斯权:

通过求解M的最小特征值的特征向量,获得p点的法向量,其中参数σ可以选择跟点云密度相关的参数。

具体σ的计算方法,我是用了paper:Efficient Simplification of Point-Sampled Surfaces,的相关方法进行计算的,代码如下:

  1. int vn=m_Tmesh->vertices.size();
  2. if(vn==m_Radius.size())return;
  3. m_Radius.resize(vn);
  4. for (int i=0;i<vn;i++)
  5. {
  6.     float max_dist=0;
  7.     for (int j=0;j<m_Tmesh->neighbors[i].size();j++)
  8.     {
  9.         float nei_dist=len2(m_Tmesh->vertices[m_Tmesh->neighbors[i][j]]-m_Tmesh->vertices[i]);
  10.         if (nei_dist>max_dist)
  11.         {
  12.             max_dist=nei_dist;
  13.         }
  14.     }
  15.     m_Radius[i]=sqrt(max_dist)/3;//估计方法见Efficient Simplification of Point-Sampled Surfaces
  16. }

上面的m_Radius代表的就是σ了。于是根据上面点云每个顶点法矢的求解公式,可以写出代码如下:

  1. for (long i = 0; i < nv; i++)
  2. {
  3.     //计算C矩阵,也就是协方差矩阵
  4.     double C[3][3] = { {0,0,0}, {0,0,0}, {0,0,0} };
  5.     for (int ni=0; ni<m_Tmesh->neighbors[i].size();ni++) //遍历点云的k近邻顶点
  6.     {
  7.         int ind =m_Tmesh->neighbors[i][ni];
  8.         vec d = vertices[ind]-vertices[i];
  9.         float weight_distance=exp(-len2(d)/(m_Radius[i]*m_Radius[i]));//m_Radius[i]是点云当前点i的密度
  10.         for (int l = 0; l < 3; l++)
  11.             for (int m = 0; m < 3; m++)
  12.             {
  13.                 C[l][m] +=weight_distance*d[l] * d[m];//高斯权重
  14.             }
  15.     }
  16.     //对协方差矩阵C进行特征值分解
  17.     Array2D<double> jz(3,3);
  18.     for (int l=0; l<3; l++)
  19.     {
  20.         for(int j=0;j<3;j++)
  21.         {
  22.             jz[l][j]=C[l][j];
  23.         }
  24.     }
  25.     Eigenvalue<double> aa(jz);//特征值分解可以调用Eigen库比较方便
  26.     Array2D<double> CC(3,3);
  27.     //获取分解结果,协方差矩阵的三个特征向量
  28.     aa.getV(CC);
  29.     m_Tmesh->normals[i]=vec(CC[0][0],CC[1][0],CC[2][0]);//第一列向量就是最小特征值对应的特征向量
  30.     if(len2(m_Tmesh->normals[i])>1e-10)normalize(m_Tmesh->normals[i]);//法矢归一化
  31. }

C、法矢方向一致化


定义成本函数:

通过上 面的特征向量的方法,求解出来的法矢,其实是不具有方向统一性的。假设每个点的最小特征值对应的特征向量,然后归一化后,法矢为V,然而其实:-V也是法 矢。因此我们需要对所以得顶点,规定一个统一的方向,比如让所有的法矢朝向曲面的外部。因此需要做统一的法矢方向归一化处理,这一步也叫法矢方向一致化, 主要用MST算法。

算法首先为点云所以顶点,定义一个成本函数:

d为从点p指向点q的单位向量,np、nq分别为点p、q的法矢,p、q是邻接顶点,这一步的代码如下:

  1. // 3a. Compute a cost between every pair of neighbors for the MST
  2. // cost = 1 – |normal1.normal2|
  3. //
  4. vector< vector<float>>costs(nv);
  5. for (int i=0;i<nv;i++)
  6. {
  7.     costs[i].resize(m_Tmesh->neighbors[i].size());
  8. }
  9. bool cost_method=true;
  10. pragma omp parallel for
  11. for(int i=0;i<nv;i++)
  12.    {
  13.     int m=m_Tmesh->neighbors[i].size();
  14.     for(int j=0;j<m; j++)
  15.     {
  16.         //文献Consolidation of Unorganized Point Clouds for Surface Reconstruction 4.1的公式
  17.         float ndot=m_Tmesh->normals[i] DOT m_Tmesh->normals[m_Tmesh->neighbors[i][j]];
  18.         float f;
  19.         if (cost_method==false)
  20.         {
  21.             vec vpq1=vertices[m_Tmesh->neighbors[i][j]] – vertices[i];
  22.             normalize(vpq1);
  23.             f=fabs(vpq1 DOT m_Tmesh->normals[i])+fabs(vpq1 DOT m_Tmesh->normals[m_Tmesh->neighbors[i][j]]);
  24.         }
  25.         else
  26.             f= 1.0- fabs(m_Tmesh->normals[i] DOT m_Tmesh->normals[m_Tmesh->neighbors[i][j]]) ;
  27.         costs[i][j]=f;
  28.     }
  29.    }

成本函数有的时候,在有的其他paper中也可以定义为:cost = 1 -|np.nq|

有了成本函数以后,接着就要以这个为法矢方向调整的依据了,进行广度优先遍历。

MST算法实现:

首先从点云模型中,选择z值最大的点,作为广度优先遍历的种子点,调整种子点的法矢方向,使得其与向量(0,0,1)的夹角小于0,这样可以保证最后所有的点云的法矢指向曲面的外侧。

接着以成本函数为权值,遍历点云广度优先遍历点云,由i点遍历到其邻接顶点j时,如果:

ni*nj<0

那么将j点的法矢调整方向,否则nj方向不需要调整。代码实现如下:

  1. // 3b. 法矢方向一致化
  2. int orientationPropagation=1;
  3. if(orientationPropagation)
  4. {
  5.     vector<int> nearby ; //未被访问的邻接顶点
  6.     int first=0; //选择索引为0的点,作为遍历的种子点
  7.     isVisited.resize(nv);
  8.     isVisited[first] = 1;
  9.     nearby.reserve(m_Tmesh->neighbors[first].size());
  10.     for(int j=0;j<m_Tmesh->neighbors[first].size();j++)
  11.     {
  12.         int nid=m_Tmesh->neighbors[first][j];
  13.         nearby.push_back(nid);
  14.     }
  15.     double cost,lowestCost;
  16.     int cheapestNearby = 0, connectedVisited = 0;
  17.     // 直到所有的点全部被遍历完毕
  18.     while(nearby.size()>0)
  19.     {
  20.         // 对于每个邻接顶点而言:
  21.         int iNearby,iNeighbor;
  22.         lowestCost = 1.0e+100;
  23.         for(int i=0; i<nearby.size(); i++)
  24.         {
  25.             iNearby = nearby[i];
  26.             for(int j=0;j<m_Tmesh->neighbors[iNearby].size();j++)
  27.             {
  28.                 iNeighbor = m_Tmesh->neighbors[iNearby][j];
  29.                 if(isVisited[iNeighbor])
  30.                 {
  31.                     cost = costs[iNearby][j];
  32.                     if(cost<lowestCost)
  33.                     {
  34.                         lowestCost = cost;
  35.                         cheapestNearby = iNearby;
  36.                         connectedVisited = iNeighbor;
  37.                         if(lowestCost<0.05)
  38.                         {
  39.                             i = nearby.size();
  40.                             break;
  41.                         }
  42.                     }
  43.                 }
  44.             }
  45.         }
  46.         //法矢点乘小于零 那么需要方向调整
  47.         if((m_Tmesh->normals[cheapestNearby] DOT m_Tmesh->normals[connectedVisited])<-1e-10)
  48.         {
  49.             m_Tmesh->normals[cheapestNearby]=(-1.0f)*m_Tmesh->normals[cheapestNearby];
  50.         }
  51.         isVisited[cheapestNearby]= 1;
  52.         //移除已然被访问的点
  53.         vector<int>::iterator iter;
  54.         iter=find(nearby.begin(),nearby.end(),cheapestNearby);
  55.         nearby.erase(iter);
  56.         // 加入未被访问的顶点
  57.         for(int j=0;j<m_Tmesh->neighbors[cheapestNearby].size();j++)
  58.         {
  59.             iNeighbor = m_Tmesh->neighbors[cheapestNearby][j];
  60.             if(isVisited[iNeighbor] == 0)
  61.             {
  62.                 vector<int>::iterator iter1;
  63.                 iter1=find(nearby.begin(),nearby.end(),iNeighbor);
  64.                 if (iter1==nearby.end())
  65.                 {
  66.                     nearby.push_back(iNeighbor);
  67.                 }
  68.             }
  69.         }
  70.     }
  71.     nearby.clear();
  72. }

当然求得点云法矢后,建议进行法矢平滑处理,这样最后重建的效果会好一些,具体可以参考我的另外一篇博文。

3、对点云的最小包围盒(bbox)进行体素网格划分。

这一步不解释,知道三维体素的人都懂,基础知识,具体体素要划分的多细得看你的需求, 可以用点云模型的bbox,进行x,y,z轴三个方向的划分,比如可以各个方向划分1000个刻度,当然刻度多大,还是得看点云模型的尺度大小的,或者你 可以先求出点云的密度的统计,在根据密度的大小进行划分。

4、求取体素的八个顶点的有向距离场

将切平面作为待重建曲面M,可以构造体素点到M的有向距离函数为:

xi为点云模型顶点,与p的最近点,ni为i点对应的法矢。这样就相当于求取p点到曲面的近似最短距离了,当然这个距离是有向的,如果点p位于M的外面,那么它的有向距离应该是正的,否则为负值。

也就是说这一步我们需要计算每个体素顶点p到曲面的最短距离,而最短的距离的计算,就是通过计算p的最近点xi,xi为点云顶点,p到xi的最短距离,就相当于其到xi切平面的距离。

  1. //计算有符号距离 
  2. void DataSets::signed_distance()
  3. {
  4.     vector<vec>& vertices=m_Tmesh->vertices;
  5.     int vn= vertices.size();
  6.     if(!vn) return;
  7.     for(long i=0;i<3;i++)
  8.     {
  9.         bounds[i*2]-=this->SampleSpacing*1.5;
  10.         bounds[i*2+1]+=this->SampleSpacing*1.5;
  11.     }
  12.     topleft[0] = bounds[0];
  13.     topleft[1] = bounds[2];
  14.     topleft[2] = bounds[4];
  15.     bottomright[0] = bounds[1];
  16.     bottomright[1] = bounds[3];
  17.     bottomright[2] = bounds[5];
  18.     for(int i=0;i<3;i++)
  19.     {
  20.         dim[i] = int((bottomright[i]-topleft[i])/this->SampleSpacing+1);
  21.     }
  22.     float   *v0 = &vertices[0][0];
  23.     KDtree  *kd = new KDtree(v0, vn);
  24.     sd=new float[dim[0]*dim[1]*dim[2]];          //存储各点的有符号距离
  25.     m_direction_field.clear();
  26.     m_direction_field.resize(dim[0]*dim[1]*dim[2]);
  27.     float   radius=16.0*SampleSpacing;
  28.     float   radius1=3.0*rou;
  29.     int count_cute=0;
  30.     for(int z=0;z<dim[2];z++)
  31.     {
  32.         for(int y=0;y<dim[1];y++)
  33.         {
  34.             for(int x=0;x<dim[0];x++)
  35.             {
  36.                 float   xx = topleft[0] + x*this->SampleSpacing;
  37.                 float   yy=topleft[1] + y*this->SampleSpacing;
  38.                 float   zz = topleft[2] + z*this->SampleSpacing;;
  39.                 long    iClosestPoint;
  40.                 vec pointp(xx,yy,zz);
  41.                 const float *match = kd->closest_to_pt(pointp, radius);
  42.                 if(match) iClosestPoint = (match – v0) / 3;
  43.                 else
  44.                 {
  45.                     sd[count_cute]=UNDEF;
  46.                     m_direction_field[count_cute]=vec(UNDEF,UNDEF,UNDEF);
  47.                     count_cute++;
  48.                     continue;
  49.                 }
  50.                 vec vectorp=pointp-vertices[iClosestPoint];
  51.                 vec plane_normal(m_Tmesh->normals[iClosestPoint][0],m_Tmesh->normals[iClosestPoint][1],m_Tmesh->normals[iClosestPoint][2]);
  52.                 float dist_plane=vectorp DOT plane_normal;
  53.                     sd[count_cute]=dist_plane;
  54.                     m_direction_field[count_cute]=dist_plane*plane_normal;
  55.                     count_cute++;
  56.             }
  57.         }
  58.     }
  59.     delete kd;
  60.  }

5、根据体素的八个顶点的符号,结合MC算法,确定求交关系。

因为我们上面求解的是有向距离场,也就是说每个体素有八个顶点,其八个顶点通过上一步我们可以求得有向距离场。那么MC算法的等值面的提取原理是什么呢?

其实如果一个体素,与待重建的曲面相交,那么该体素的8个顶点的有向距离,肯定会有正有负。因此我们只需要处理那种8个顶点中有正有负的体素就可以了。而MC算法就是根据相交情况,进行查表,得到相交的三角网格模型的。

首先对体素的8个顶点进行分类,以判断其顶点是位于等值面之外,还是位于等值面之内。再根据8个顶点的状态,确定等值面的剖分模式。顶点分类规则为:

1)如体素顶点的数据值大于或等于等值面的值,则定义该顶点位于等值面之外, 记为正点,即“1“

2)如体素顶点的数据值小于等值面的值,则定义该顶点位于等值面之内,记为负点, 即“0″

由 于每个体素共有8个顶点,且每个顶点有正负两种状态,所以等值面可能以 =256种方式与一个体素相交。通过列举出这256种情况,就能创建一张表格,利用它可以查出任意体素中的等值面的三角面片表示。如果考虑互补对称性,将 等值面的值和8个角点的函数值的大小关系颠倒过来,即将体素的顶点标记置反(0变为1, 1变为0),这样做不会影响该体素的8个角点和等值面之间的拓扑结构,可将256种方式简化成128种。其次,再利用旋转对称性,可将这128种构型进一 步简化成15种。图3.2给出了这15种基本构型[131其中黑点标记为“1”的角点。

根据8个顶点有向距离场,的正负情况进行查表,可以得到三角面片。然后三角面片的顶点的坐标,也就是体素点边与待重建曲面的相交点,这个可以用体素点边的两个顶点进行线性插值得到,网格曲面每个顶点的坐标。

于是有了网格曲面的每个顶点,也有了三角面片,这样网格曲面也就重建完毕了。OK算法至此结束。

参考文献:

1、《基于点云的鞋楦三角网格曲面重构》

2、《点云模型法矢优化调整》

转载注明来源:CV视觉网 » 图形处理(九)点云重建(下)法矢求取、有向距离场等值面提取

分享到:更多 ()
扫描二维码,给作者 打赏
pay_weixinpay_weixin

请选择你看完该文章的感受:

0不错 1超赞 0无聊 0扯淡 0不解 0路过

评论 3

评论前必须登录!