[發明專利]一種基于區域配準的深度點云三維重建方法

在審

• 申請號: 201611233141.3
• 申請日: 2016-12-28
• 公開（公告）號: CN106651926A
• 公開（公告）日: 2017-05-10
• 發明（設計）人: 全紅艷;趙含放
• 申請（專利權）人: 華東師范大學
• 主分類號: G06T7/33
• 分類號: G06T7/33;G06T7/593
• 代理公司: 上海藍迪專利商標事務所(普通合伙)31215
• 代理人: 徐筱梅,張翔

1.一種基于區域配準的深度點云三維重建方法，其特征在于，該方法包括以下具體步驟：步驟1：對彩色圖像進行區域分割

將彩色圖像每個像素初始化為圖G的頂點，像素之間的顏色差異初始化為頂點之間的邊權值；計算所有像素之間的顏色差異，計為圖G的邊權值h，按照公式(1)計算：

$h_{ij}=\sqrt{{(r_i-r_j)}^2+{(g_i-g_j)}^2+{(b_i-b_j)}^2}---\left(1\right)$

其中，h_ij是像素點p_i與p_j之間的邊權值，(r_i,g_i,b_i)與(r_j,g_j,b_j)分別是兩個像素點的紅綠藍三個顏色值；為了進一步進行區域分割，先將圖G中每一個頂點當作一個區域，然后根據公式(2)對每兩個區域進行聚類分析，如果兩個區域滿足公式(2)，那么將他們合并為一塊新區域；

Dif(C_i，C_j)≤min(Int(C_i)+f/|C_i|,Int(C_j)+f/|C_j|) (2)

其中，C_i和C_j是待分割的兩個區域，Dif(C_i，C_j)是連接兩個區域C_i和C_j之間的最小的邊權值，Int(C_i)和Int(C_j)分別為區域C_i和C_j包含的圖所對應的最小生成樹中最大的邊權值；|C_i|和|C_j|分別是區域C_i和C_j中的像素點個數,f是分割的閾值參數，設為50；對所有像素點進行處理后，彩色圖像被過分割成多個區域；

步驟2：對相似的過分割區域進行合并

首先，根據公式(3)計算每個區域的協方差矩陣A(u(S))：

$A\left(u\right(S\left)\right)=\frac{1}{\left|u\left(S\right)\right|}{\underset{x\∈u(S}{\Σ}}_{)}(x-c(S\left)\right){(x-c(S\left)\right)}^T---\left(3\right)$

其中，c(S)是區域S的中心點的三維坐標，u(S)是區域S對應的三維點云；|u(S)|為u(S)的頂點數量；x是點云上u(S)任意一個點的三維坐標；對A(u(S))進行矩陣的特征分解，得到的三個特征值從大到小分別為w₁(S)、w₂(S)、w_n(S)，對應的三個特征向量分別為T₁(S)、T₂(S)、N(S)；

然后，對于圖像中任意兩個相鄰區域S_i,S_j,按照公式(4)計算它們之間的相似度c_smooth：

$c_{smooth}=\frac{\max \left\{\right||N\left(S_i\right)\·\left(c\right(S_i)-c(S_j\left)\right)||,||N\left(S_j\right)\·\left(c\right(S_i)-c(S_j))\left|\right|\}}{\left|\right|N\left(S_i\right)\·N\left(S_j\right)\left|\right|}---\left(4\right)$

其中||·||表示數量積計算，N(S_i)與N(S_j)是兩個區域計算得到的最小的特征向量；如果兩個區域的相似度c_smooth小于閾值常數c_thod，相鄰區域S_i、S_j將合并成為新區域；其中c_thod設為0.02m；對相鄰時刻所有區域都進行相似度計算和合并后，得到相鄰時刻兩組合并后的區域集合；

步驟3：匹配相鄰時刻的兩組區域，對所有匹配的區域所對應的點云進行三維配準

從時刻t_k的區域集合中取一塊區域S_m，再從下一時刻t_k+1的區域集合中取一塊區域S_n，組成區域對(S_m,S_n)；然后按照公式(5)，計算S_m與S_n之間的空間位置關系r_mn；

$r_{mn}=w_{mn}{d}_{mn}\frac{\left|\right|N\left(S_m\right)\left|\right|\·\left|\right|N\left(S_n\right)\left|\right|}{\left|\right|N\left(S_m\right)\·N\left(S_n\right)\left|\right|}---\left(5\right)$

其中，d_mn為兩個點云區域之間所有最臨近點對的平均歐氏距離，按照公式(6)計算；

$d_{mn}=\frac{1}{\min \left\{\right|u\left(S_m\right)|,|u\left(S_n\right)\left|\right\}}\underset{(x_n,y_n,z_n)\∈u\left(S_n\right)}{\underset{(x_m,y_m,z_m)\∈u\left(S_m\right)}{\Σ}}\sqrt{{(x_m-x_n)}^2+{(y_m-y_n)}^2+{(z_m-z_n)}^2}---\left(6\right)$

其中，(x_m,y_m,z_m)是u(S_m)上的任意點的三維坐標，(x_n,y_n,z_n)是在u(S_n)上距離(x_m,y_m,z_m)歐式距離最小的點的三維坐標；

w_mn表示兩個點云的重疊率，按照公式(7)進行計算；

$w_{mn}=\frac{N_w}{\min \left\{\right|u\left(S_m\right)|,|u\left(S_n\right)\left|\right\}}---\left(7\right)$

其中，N_w是兩個點云中歐式距離最小的點對距離小于d_thod的點對數，其中d_thod設為0.01m；

然后，對相鄰兩個時刻所有區域的點云進行權值初始化，具體地，區域S_m和S_n的權值e(S_m)，e(S_n)分別初始化為S_m和S_n對應點云u(S_m)和u(S_n)的頂點個數占該時刻所有點云的頂點個數合的比例；如果r_mn小于閾值r_thod，那么將該組點云區域對(S_m,S_n)加入預處理集合M中，并分別按照公式(8)更新權值e(S_m)，e(S_n),其中r_thod設為0.02；

$\begin{array}{c}e\left(S_m\right)\←e\left(S_m\right)-\min \{e\left(S_m\right),e\left(S_n\right)\}\\ e\left(S_n\right)\←e\left(S_n\right)-\min \{e\left(S_m\right),e\left(S_n\right)\}\end{array}---\left(8\right)$

對所有區域對的點云進行權值計算，迭代更新所有點云的權值，最終得到包含所有區域對的預處理集合M；

進一步對預處理的區域對進行對應點云的幾何三維配準；首先從M中取出一對區域，對這兩塊區域對應的三維點云分別進行幾何特征提取；具體地，對點云u(S_m)和u(S_n)中的每個點分別計算FPFH特；對點云每個點周圍距離小于R的點進行特征統計形成點特征直方圖，其中特征統計包括對距離、角度的統計，R取值為0.5m；然后利用采樣一致性初始配準算法計算物體的四階變換矩陣，具體地，先在u(S_m)中尋找采樣點，然后在另一個點云u(S_n)中查找與采樣點FPFH特征相似的所有點，計算每個點對的四階變換矩陣，選擇其中誤差最小的剛體變換矩陣T_kp作為最終變換矩陣；k表示第k時刻，p表示在預處理集合M中的第p對匹配區域；其中，采樣一致性初始配準算法中的采樣點之間的距離不應小于閾值d_min，設為0.1m；

步驟4：融合通過三維配準得到的所有區域的變換矩陣

對于同一個時刻k，將該時刻所有區域的變換矩陣T_kp融合成該時刻的完整變換矩陣T_k，假設時刻k總共有p_k組區域對，按照公式(9)計算融合后的變換矩陣T_k；

$T_k=X\·\left(\begin{array}{c}\frac{1}{\underset{p=1}{\overset{p_k}{\Σ}}{w}_{kp}}\underset{p=1}{\overset{p_k}{\Σ}}{w}_{kp}{T}_{kp}\\ \frac{1}{\underset{p=1}{\overset{p_k}{\Σ}}{d}_{kp}}\underset{p=1}{\overset{p_k}{\Σ}}{d}_{kp}{T}_{kp}\end{array}\right)---\left(9\right)$

其中X為權值矩陣，設為(0.70.3)，T_kp是時刻k時的第p對區域點云配準得到的變換矩陣，w_kp是時刻k時的第p對區域點云的重疊率，按照公式(7)進行計算，d_kp是時刻k時的第p對區域點云的平均歐式距離，按照公式(8)進行計算；

步驟5：將不同時刻的深度點云進行點云融合與三維表面重建

將空間劃分成128³個柵格立方體，將點云分布在柵格立方體中，然后計算每幀融合后每個柵格立方體的TSDF值；

為了計算每個柵格立方體中的截斷符號距離函數，首先要計算每個柵格立方體中的符號距離函數sdf_kg(x)，k表示第k時刻，g表示第g個柵格立方體，sdf_kg(x)代表著點云上任意一點x的符號距離函數，按照公式(10)進行計算：

sdf_kg(x)＝dis(v_g,t_k)-D_k (10)

其中，t_k是相機光心在三維坐標系下的坐標，相機光心坐標設為(320,240)，v_g是柵格g在三維坐標系下的坐標，dis(v_g,t_k)表示柵格到光心的距離，D_k是圖像中x對應的像素點所對應的深度值，由Kinect設備相機采集得到；如果sdf_kg(x)大于0，tsdf_kg(x)按照公式(11)計算，否則按照公式(12)計算：

tsdf_kg(x)＝min(1,sdf_kg(x)/m_truncation)(11)

tsdf_kg(x)＝min(-1,sdf_kg(x)/m_truncation)(12)

其中m_truncation為常量，取值為10；然后計算權重wt_kg(x)，初始值為1；W_kg(x)為迭代后的權重，計算公式如(13)所示；

W_kg(x)＝W_k-1,g(x)+wt_kg(x) (13)

最后記錄迭代后的加權TSDF值，標記為TSDF_kg(x)，按照公式(14)計算；

${\mathrm{TSDF}}_{kg}\left(x\right)=\frac{{\mathrm{TSDF}}_{k-1,g}\left(x\right)W_{k-1,g}\left(x\right)+{\mathrm{tsdf}}_{kg}\left(x\right)w_{kg}\left(x\right)}{W_{k-1,g}\left(x\right)+w_{kg}\left(x\right)}---\left(14\right)$

將W_kg(x)和TSDF_kg(x)存儲在對應的柵格g中，進行下個柵格g+1的計算；當對所有時刻的柵格都計算了截斷符號距離函數后，所有點云已被存儲在柵格立方體中；然后使用移動立方體算法抽取等值面，繪制三維表面，得到立體模型。