基于圖像的三維重建典型技術綜述

孫海迅，羅健欣，張艷艷，鄭義桀，潘志松

（中國人民解放軍陸軍工程大學 指揮控制工程學院，江蘇 南京 210001）

0 引言

三維重建是指在計算機中建立現實物體或場景三維模型的關鍵技術。隨著計算機視覺的發展，三維重建在醫學影像判讀、文物修復、自動駕駛、虛擬現實等領域有著廣泛應用。按照重建方法，可將其分為傳統方法和深度學習方法兩大類。傳統方法主要是對于在不同視角獲取的物體照片，利用相關算法，恢復其三維結構。深度學習方法是通過訓練神經網絡預測物體的三維模型［1］。深度學習方法對算力要求較高，而傳統基于圖像的三維重建價格低廉且算法相對成熟，因而應用更加廣泛。

基于圖像的三維重建主要包括以下流程：圖像特征點檢測與匹配、稀疏點云重建、稠密點云重建、三維表面重建、表面紋理貼圖。本文針對以上流程中的關鍵環節，選取幾種代表性方法進行綜述。

1 特征點檢測與匹配

特征點是圖像中亮度變化劇烈的點或圖像邊緣曲線上曲率的極大值點，它同周圍的鄰近點有著明顯差異［2］，目前已有多種特征點檢測方法［3］。在三維重建過程中，可以通過對不同視角下圖像特征點的檢測與匹配，求解相機姿態，進而獲取三維點云。本文重點介紹兩種經典方法——Harris角點檢測［4］和SIFT 特征點檢測［5］。

1.1 Harris角點檢測

角點指圖像上兩個方向上邊緣線的接合點。在一幅圖像上取一小塊特征區域，如果區域的灰度沿各方向都有較大變化，則認為窗口內存在角點。

設（x，y）為圖像上某處像素坐標，I（x，y）代表該點處對應的光度，（u，v）代表位移，I（x+u，y+v）代表位移后的光度，w（x，y）代表權重函數，本文選擇高斯函數，則可用E（u，v）度量特征區域位移后與位移前的光度差異，以此判斷區域內是否包含角點。

其中，對I(x+u，y+v)進行一階泰勒展開，即：

其中，Ix、Iy表示圖像在(x，y)坐標點處的亮度分別對x、y的偏導數，因而可將上式重寫為：

通過對M矩陣的特征值進行分析即可判斷Ω 區域是否包含角點，當M的兩個特征值λ1、λ2都比較大時，則認為Ω 區域包含角點。同時，文獻［4］給出了如下判別方法，當R大于某閾值，認為Ω 區域存在角點。

1.2 SIFT特征點檢測

文獻［5］提出了一種對旋轉、縮放、明暗保持不變的特征點檢測方法，通過高斯差分算子（Difference of Gaussian，簡稱DOG 算子）尋找特征點。高斯空間與高斯差分空間如圖1所示。

Fig.1 Gaussian space and Gaussian difference space圖1 高斯空間與高斯差分空間

首先對原始圖像進行降采樣形成不同的組（octave），在每組內對圖像用不同的高斯核函數（標準差分別為σ、kσ、k2σ......）做卷積形成多張圖像。計算方法為：

其中，I(x，y)表示原圖像，G(x，y，σ)表示標準差為σ的核函數，L(x，y，σ)表示卷積之后生成的一幅高斯圖像，以此構建高斯圖像空間（Gussian Space）。

同一組內相鄰兩幅圖像相減得到高斯差分圖像空間（DoG Space），計算方法為：

D(x，y，σ)表示生成的一幅高斯差分圖像，并且每一組內有效差分數S與高斯圖像空間的層數N滿足如下關系：

在高斯差分圖像空間求極值得到圖像特征點，即對每個點與它所在圖像與上下相鄰兩層圖像周圍26 個像素點的亮度進行比較，如果大于或者小于這26 個點，則作為極值點。

進一步地，通過在亞像素位置上尋找極值點的更精確位置能夠從本質上改善特征點的匹配穩定性［5］，圖像極值點如圖2 所示。于是，在x0=(x0，y0，σ0)T處對函數D(x，y，σ)進行泰勒展開得：

Fig.2 Extremum point of the image圖2 圖像極值點

令δx=x-x0，將D(x)對其求偏導使之為0，得：

于是精確的極值位置為x=x0+δx，對應的函數值：

其中，偏導數矩陣?DT(x0)、?2D(x0)是用像素點附近點的差分形式近似表示［5］。

此外，類似于Harris 角點檢測的方法，對圖像平面上x處的二階Hessian 矩陣的特征值施加一定約束，能剔除某些只在一個方向上與其他點有明顯差異的點，使得檢測到的特征點更加穩定。約束條件為：

圖像的主方向通過統計梯度直方圖的方法確定，將圖像的x軸旋轉到主方向上，算法才具有旋轉不變性。

SIFT 特征描述子（128 維）同樣是統計局部區域梯度信息而生成，通過最近鄰搜索方式進行匹配，并將最近鄰與次近鄰的距離之比限定在較小范圍內，以此得到較為精確的匹配對。

SIFT 特征點檢測比Harris 角點檢測更魯棒，對旋轉、尺度縮放、亮度變化能保持不變，但實時性不高、運算速度較慢，因而出現了諸如SURF［6-7］等方法的改進。

2 稀疏點云重建

稀疏點云重建是根據圖像上特征點的匹配關系，通過運動恢復結構（SFM）生成稀疏點云的方法。運動恢復結構是指通過圖像間的對應關系，重構相機外部參數（位置和姿態）和內部參數（焦距和徑向畸變）。在實現SFM 的過程中，通常包含三角測量，以此求解三維點坐標，得到稀疏點云。

2.1 基本概念

2.1.1 針孔相機模型

針孔相機模型的坐標系包括：世界坐標系、相機坐標系、歸一化像平面坐標系和物理像平面坐標系。世界坐標系是物體所在現實世界的坐標系，實際計算時常用相機的初始位置建立世界坐標系；相機坐標系是在以相機中心為坐標原點建立的坐標系，其Z軸通常指向要拍攝的物體。世界坐標系到相機坐標系的坐標變換如圖3 所示，相機坐標系下的坐標變換如圖4所示。

Fig.3 Coordinate transformation from world coordinate system to camera coordinate system圖3 世界坐標系到相機坐標系的坐標變換

Fig.4 Coordinate transformation in camera coordinate system圖4 相機坐標系下的坐標變換

設世界坐標系下某一點Xw（xw，yw，zw），在相機坐標系下坐標為Xc（xc，yc，zc），二者可通過旋轉平移變換得到，于是滿足以下關系：

寫成齊次形式：

歸一化像平面是人為設置的虛擬平面，將其z坐標（距離相機光心距離）設為1。設點Xc的歸一化坐標為（x，y），則：

物理像平面平行于歸一化平面，指相機CCD 陣列所在的平面。物理像平面坐標也即圖像的像素坐標，歸一化坐標乘上焦距f和分辨率α（可以合并寫為fα），再進行坐標平移后（因為物理像坐標通常以圖像左上方為原點）即可得到物理像平面坐標。

因此，從空間點的世界坐標Xw到圖像坐標的變換關系為：

其中，K為內參矩陣，［R t］為外參矩陣。

2.1.2 徑向畸變矯正

相機的針孔模型與真實的相機成像不盡相同，真實的相機會由于透鏡形狀引起圖像畸變——徑向畸變。因此，要從二維圖像中準確提取度量信息，必須進行攝像機標定。Zhang［8］提出一種利用棋盤格標定相機的方法，標定過程僅需打印棋盤格，并從不同方向拍攝圖片，不僅靈活方便，而且魯棒性好，如圖5所示。

Fig.5 Comparison before and after correction of radial distortion圖5 徑向畸變矯正前后對比

該方法假設畸變后的歸一化坐標與畸變之前坐標存在如下關系：

其中，k1、k2為徑向畸變系數，r為坐標點到歸一化平面中心點的半徑。畸變后的圖像坐標為：

在理想情況下：

聯立以上等式，可得：

其中，、是通過觀察棋盤格上已知的某點在未矯正圖像上的位置而得到，u、v是理想情況下利用前面的投影關系計算得到，r是該點到棋盤中心的距離。在棋盤格上選定多點，即可利用最小二乘法解方程得到k1、k2。之后，在矯正后的圖像上，即可利用如下關系得到點（u，v）對應的、，并將矯正前坐標（、）處的顏色值作為矯正后圖像上（u，v）坐標的顏色值，生成矯正后的圖像。

2.2 相機姿態恢復

在對極幾何［9］中，常通過圖像坐標點的對應關系，求解基本矩陣或單應性矩陣進而求解相機姿態。

2.2.1 通過基本矩陣與本質矩陣求解相機姿態

如圖6 所示，X為一空間三維點，O1、O2為兩個相機中心，則XO1O2所在的平面稱為極平面，π1、π2為兩個相機的物理像平面，x1、x2分別為X在兩個物理像平面的像素坐標。設、分別為x1、x2對應的歸一化平面坐標，[R，t]為O2坐標系相對于O1坐標系的坐標變換矩陣，兩個相機的內參矩陣分別為K1、K2，則有如下約束關系：

Fig.6 Diagram of epipolar geometry圖6 對極幾何示意圖

在圖像特征點檢測與匹配的基礎上，利用以上約束關系，即可求得基本矩陣和本質矩陣，進而得到相機的外參數（相機姿態）。文獻［10］給出了利用8 對在圖像上匹配的特征點計算基本矩陣F的方法，稱為8點法。具體為：

設x1=(u1，v1，1)T，x2=(u2，v2，1)T為一對匹配的特征點，根據約束1=0，即：

若令f=(F11，F12，F13，F21，F22，F23，F31，F32，F33)T，則有：

也即每一對匹配點能對F矩陣提供一個約束，而若令F的任一元素歸一化（所有元素同時除以某一元素），并不影響1=0 的成立，即F有8 個自由度，因此至少需要8對匹配的特征點即可求得F，當匹配的特征點大于8 時，也可得利用最小二乘法求解。

為增加8 點法的魯棒性，文獻［11］采用了將圖像坐標歸一化的方法，也稱為歸一化8 點法。文獻［12］提出要對F矩陣增加奇異值約束。文獻［13］提出了一種根據奇異值約束對矩陣F進行重構的方法，即在利用8 點法得到F的基礎上，將F進行SVD分解：

利用奇異值約束將基本矩陣重構為：F，=Udiag(σ1，σ2，σ3)VT，這樣可以使得F-F，的Frobenius 范數最小。

在實際應用中，特征點可能因噪聲產生誤匹配而使得基本矩陣的估計不準確。為解決這一問題，文獻［14］采用隨機抽樣一致性（RANSAC［15］）方法估計基本矩陣，具體為：

（1）計算采樣次數。

其中，z代表期望采樣達到的成功率，M代表為達到這一成功率至少需要的采樣次數，k代表每次采樣求解基本矩陣所需要的點數，p代表所有樣本點中內點（匹配準確的特征點）的概率。

（2）對匹配的特征點集中隨機采樣8 對，利用8 點法估計初始矩陣F。

（3）用Sampson 距離［16］度量點集中所有點，將點劃分為內點和外點（誤匹配的特征點），并記錄當前計算所得F矩陣對應的內點個數。一對特征點的Sampson 距離可表示為：

其中，x1、x2表示不同圖像上的一對特征點，F為步驟（2）計算得到的矩陣，下標x、y分別表示某一向量的x分量和y分量。

（4）重復步驟（2）（3）M次，選擇內點個數最多的F矩陣和相應的內點。

（5）對所有內點重新執行步驟（2）得到基本矩陣F’。

相機姿態的恢復需要求解本質矩陣。本質矩陣E是一個秩為2 的矩陣，具有5 個自由度，有兩個相等的奇異值和一個0 奇異值［17］。與求解基本矩陣F所用的8 點法類似，本質矩陣的求解可以采用5 點法［18］，也可在求得基本矩陣F的基礎上，利用以下關系求本質矩陣初始解：

同樣地，要對E施加奇異值約束并重構［9］：

在文獻［9］中證明，SVD 分解不唯一，如果可以分解為式（30），則也一定可以分解為：E=Udiag(1，1，0)VT。如果可以將本質矩陣進行SVD 分解，則E=[t]×R有如下4 種可能情況：

從圖7 可以看出，4 種情況中只有一種情況是空間點X同時位于兩個相機的鏡頭前方，因而只要用一個點做測試，看它是否位于兩個相機的前方，即可選擇出正確的相機姿態。

Fig.7 Choice of camera pose圖7 相機姿態選擇

2.2.2 用單應性矩陣求相機姿態

文獻［19］指出，當圖像上匹配的特征點均來自于空間中同一個平面時，無論有多少這樣的匹配點，它們對基本矩陣F的約束不會多于4 個，這將導致求解基本矩陣F的結果不穩定。此外，當相機姿態只有旋轉沒有平移時（從O1變換到O2是純旋轉，t=0），E=[t]×R=0，難以通過對E的分解求解R。

如圖8 所示，若空間中一平面π 方程為：nTX+d=0，則空間點X 對應的兩幅圖像上的投影x1、x2有如下關系：

Fig.8 Diagram of homography transformation圖8 單應性變換示意圖

其中K1、K2分別為兩個相機的內參矩陣。若相機為純旋轉變換，即t=0，則H=K2RK1-1。對于H的求解方法，同樣可以采用類似于求解基本矩陣F的方法，利用兩幅視圖上匹配點的約束關系x2=Hx1，H矩陣有8 個自由度，而一對匹配點能提供兩個約束，因此至少要4 對匹配點即可求解H［20］。此外，可以采用RANSAC 的方法估計H，然后便可對H進行分解以求解R和t［20］。

2.3 三角測量

在已知相機參數和匹配點的前提下，即可利用三角測量恢復空間點的三維坐標，具體方法如下：

2.3.1 直接線性變換法

文獻［9］提出了直接線性變換法。設第i個相機的投影矩陣為：Pi=Ki[Ri，ti]，空間點的三維坐標為X=[x，y，z，1]T，對應在第i個視角的圖像坐標為xi=[ui，vi，1]T，di表示X點在第i個視角下的深度，即相機坐標系下的坐標zc。則：

兩邊同時叉乘xi，得：

其中，第3 個方程與前兩個線性相關，則一個觀測視角能對X提供兩個約束，而X有3個自由度，因此要求解X，至少需要圖像上一對匹配點，也即至少要從兩個視角對其進行觀測。當點數N 大于2 時，還可以構建方程組利用最小二乘法求解：

2.3.2 RANSAC方法

在存在外點的情況下，可以用RANSAC 方法求解［15］。主要步驟如下：①RANSAC 采樣次數，設置內點閾值（重投影誤差）；②隨機采樣一對視角，計算空間點三維坐票；③計算每個視角中的重投影誤差，統計內點個數；④重復步驟②、③直到滿足采樣次數，選擇內點數最多的視角；⑤利用所有內點重新計算三維點坐標。

此外，在利用相機匹配對進行三角測量計算三維點時，可以從圖像中提取EXIF 初始值以獲取焦距［21-22］。

2.4 PnP問題

PnP 問題指根據已知的三維點坐標和對應的二維點，求解相機內外參數的問題，這在增量式SFM 中經常被用到［21］。

常用的方法為直接線性變換法，設相機的投影矩陣為：P=K[R，t]，空間點的三維坐標為X=[x，y，z，1]T，對應的圖像坐標為x=[u，v，1]T，d表示X點的深度，即相機坐標系下的坐標zc。則：

P矩陣共有12 個參數，而一對3D-2D 對應點能提供2個約束，因此至少需要6 對匹配點即可求解P。因為P[：，1：3]=KR，且K為上三角陣，而R為正交矩陣，可以對P[：，1：3]-1進行QR 分解，即可求得R-1和K-1，進而求得K和R。又P[：，4]=Kt，則t=K-1P[：，4]，得平移向量t。

文獻［23］提出一種計算PnP 問題的Kneip 算法，利用三維點構建一個平面作為過渡坐標系，進而求解相機的內外參數。此外，同樣可以用RANSAC 方法實現Kneip 算法。求解PnP 問題還有其他幾種方法，如用4 對不共面的點進行求解的P3P 方法［24］，利用4 對以上不共面或者3 對共面點進行求解的ePnP 法［25］，以及利用深度學習與傳統方法相結合的EPro-PnP 方法［26］，實現端到端的三維位姿預測。

2.5 捆綁調整

由于噪聲的影響，以上初步求得相機內外參數和三維點坐標還不夠準確，捆綁調整（Bundle Adjustment）［27］可以同時對三維點坐標和相機參數進行非線性優化。捆綁調整如圖9所示。

Fig.9 Diagram of bundle adjustment圖9 捆綁調整示意圖

設χij=1 表示第i個點在第j個相機中可見，χij=0 表示第i個點在第j個相機中不可見；Xi=(xw，yw，zw，)T表示第i個三維點的世界坐標；Pj表示第j個相機的內參以及世界坐標系與第j個相機的外參數；表示第i個點在第j個相機中的觀測點；uij表示第i個點在第j個相機中的投影點。

則在有n個三維點和m個視角時，可構建損失函數：

其中，θ=(P1，...Pm，X1，...Xn)為待優化的量。對無約束非線性最小優化問題對于最優化問題g(θ)=通常可用以下幾種方法。

2.5.1 最速下降法

假設g(θ)在θt處可微，則它在θt處有泰勒展開（展開到一階）。

其中，θ=θt+δθ，當(?g(θt))Tδθ＜0 時可保證g(θ)的值在下降，當δθ取梯度的反方向時，可達到最快下降速度。

算法流程如下：

Step1：給定初始點θ0，迭代終止閾值ε和步長λ，令t=0。

Step2：計算?g(θt)，若‖ ‖?g(θt) ≤ε停止迭代，輸出θt，否則轉下一步。

Step3：取θt+1=θt-λ?g(θt)，t=t+1，轉Step2。

2.5.2 牛頓法

最速下降法是收斂的，但最終的收斂是線性的，而且通常很慢，并且在某些情況下無法找到二階多項式的最小值，此時可用牛頓法［28］。

假設g(θ)在θt處二階可微，且假定二階導數?2g(θ)總是正定，則以它在θt處二階近似函數Q(θ)極小值作為下一次迭代點θt+1。同樣，將Q(θ)泰勒展開為：

對θ求梯度令其為0，可得：

算法流程如下：

Step1：給定初始點θ0，迭代終止閾值ε和步長λ，令t=0。

Step2：計算?g(θt)，若‖ ?g(θt) ‖≤ε，停止迭代，輸出θt，否則轉下一步。

Step3：取θt+1=θt-λ(?2g(θt))-1?g(θt)，t=t+1，轉Ste p2。

2.5.3 高斯—牛頓法

牛頓法雖然速度快但計算量大，需要保存二階Hessian 矩陣的逆矩陣，高斯—牛頓法與牛頓法類似，但區別是它基于f(θ)在θ附近的線性近似［29］，即設p(θ)=f(θ)-x，對其泰勒展開為：

令JT(θt)=?f(θt)=?p(θt)，則：

與牛頓法相比，用JT(θt)J(θt)代替?2g(θt)，降低了運算復雜度。

2.5.4 Levenberg-Marquardt法

最速下降法是用平面在局部擬合損失函數，牛頓法是用二次曲面擬合，當?2g(θ)在某區域負定時，牛頓法可能收斂到局部極大值，高斯—牛頓法在J(θ)不滿秩時可能失效［28］。Levenberg-Marquardt 法［30-31］是將最速下降法和高斯—牛頓法有機結合，使用“信賴閾”控制下降速度的方法。當收斂速度較快時，信賴閾增大，算法趨向于高斯—牛頓法；當收斂速度慢時，信賴閾減小，算法趨向于最速下降法，優點是只用到一階雅各比矩陣，計算速度快，對初始值有一定的魯棒性。

通過求解增量正規方程得到δθ。可見，當λ趨向無窮大時，(JT(θt)J(θt))δθ=-?g(θt)，近似于高斯—牛頓法；當λ趨近于0時，δθ=-?g(θt)，近似于最速下降法。

Levenberg-Marquardt 法的主要流程為：

Step4：通過求解增量正規方程，得到δθ。

2.6 運動恢復結構分類

典型的運動恢復結構方法有增量式（Incremental SFM）［32-35］、全局式（Global SFM）［36-37］、分層式［38-39］。

2.6.1 增量式運動恢復結構

主要步驟包括：圖像特征點檢測與匹配、構建圖像連接圖、初始化、增量式重建新的視角、全局捆綁調整。

（1）圖像特征點檢測與匹配。采用上文所述方法，SIFT 方法檢測的特征點比較穩定，其應用較廣泛［33-34］。

（2）圖像連接圖構建。每一張圖像形成圖上的一個結點，具有特征匹配點的兩幅圖像之間有邊相連，結點的邊越多，表示與其具有匹配的特征點的其他圖像越多。從圖10 可以看出，白天和夜晚拍攝的圖像因為光度差異會在連接圖上形成兩個簇［33］。

（3）初始化。選擇一對初始相機進行Tracks 重建（三角測量）恢復三維點坐標，然后進行Tracks 濾波和捆綁調整以優化估計的三維點坐標和相機內外參數。

其中，初始相機對的選取要遵循一定的原則：匹配點足夠多、基線足夠長、滿足單應性的匹配點盡量少［33-34］。文獻［34］指出，由于冗余數據的存在，從圖像連接圖中的密集位置進行初始化通常會使得重建更加魯棒和準確。

Tracks 濾波是指濾除Tracks 重建中生成的錯誤點（包括無窮遠處的點和重投影誤差過大的點）［21，33］，以降低后續運算復雜度。

（4）新的視角增量式重建。選擇能看到已被估計的三維點最多的視角，利用PnP 方法［21］或者RANSAC 方法［33］求解新視角的相機姿態，然后對單個相機姿態進行捆綁調整。下一步對新視角能看到的點進行Tracks 重建和濾波，以此優化三維點的坐標。

（5）全局捆綁調整。對所有已加入的相機內外參數和三維點坐標進行一次全局的非線性優化，由于這一方法運行時間較長，可以在加入多個視角之后運行一次［21］。

增量式運動恢復結構如圖11 所示。其主要優點是：對誤匹配的特征點魯棒、場景結構可在調整中不斷優化。缺點是：對初始視角的選取和視角添加順序敏感、大場景下的累計誤差會導致場景漂移、捆綁調整重復進行導致效率低。文獻［35］提出一種與RANSAC 相結合提升精度的增量式運動恢復結構方法。

Fig.11 Incremental structure from motion圖11 增量式運動恢復結構

2.6.2 全局式運動恢復結構

由于捆綁調整比較耗時，全局式運動恢復結構同時估計所有相機的旋轉矩陣和平移向量，生成三維點后，只運用一次捆綁調整進行優化。文獻［36］提出一種估計全局相機姿態的方法。

首先，在生成圖像連接圖之后，利用任意兩幅圖像的特征點匹配關系，計算兩個相機的相對旋轉矩陣Rij和平移向量tij，通過最小化如下關系求解所有相機相對于世界坐標系的旋轉矩陣。

其次，進行全局旋轉矩陣濾波，根據求得的Ri與Rj優化圖像連接圖，即當時，認為相機i與j之間存在錯誤的連接邊，要刪除。

最后，利用求得的相機姿態求解三維點并進行一次全局的捆綁調整。

文獻［37］提出一種多種信息融合的方法，在求解過程加入了GPS等輔助信息，對結果進行優化。

全局式運動恢復結構如圖12 所示。其主要優點是：沒有誤差積累，誤差均勻分布在連接圖上；對初始相機選取和相機添加順序不敏感；捆綁調整僅執行一次，重建效率高。不足之處主要體現在：相機位置求解時對匹配外點敏感，魯棒性不足；連接圖邊界被過濾，容易造成部分圖像丟失，導致重建不完整。

Fig.12 Global structure from motion圖12 全局式運動恢復結構

2.6.3 分層式運動恢復結構

文獻［38］提出基于聚類方式的分層式運動恢復結構，如圖13 所示。其主要步驟為：利用已知視角創建一個局部模型進行捆綁調整，加入一個新的視角，兩個模型融合捆綁調整。這種方法計算效率高，且算法穩定，不依賴初始相機的選取，避免了大場景中的誤差積累和漂移，但規則設計較為復雜。

Fig.13 Hierarchical structure from motion圖13 分層式運動恢復結構

3 稠密點云重建

基于圖像的三維重建中，稠密點云的獲取通常基于多視圖立體技術（Multi-view Stereo，MVS）。該技術的優點是成本低、精度較高、圖像來源廣；缺點是計算速度慢。

MVS 中常用的一個基本假設是光度一致性假設（Photo-Consistency），是指同一空間的點在不同視角的投影，光度應該相同，重建的關鍵是恢復光度一致性的點。光度一致性度量方式有：誤差平方和（SSD，Sum of Squared Difference）、歸一化交叉協方差（NCC，Normalized Cross Correlation）。具體為：

其中，f和g分別為特征描述子的特征向量、為描述子的均值向量，δf、δg為描述子向量標準差。

基于多視圖立體技術（MVS）的稠密點云重建主要有基于體素的方法、基于深度圖融合的方法和基于空間patch 擴散的方法。

3.1 基于體素的方法

基于體素的方法是指將空間劃分為很多體素，通過圖像序列的約束計算體素內三維點的方法。空間劃分包括兩種類型：一種是規則體素的劃分［40］，即將空間劃分為小立方體；另一種是不規則體素的劃分［41-42］，即將空間劃分為小四面體，如圖14所示。

Fig.14 Dense reconstruction based on voxels圖14 基于體素的稠密重建

這種方法的優點是：生成規則的點云、便于提取物體的平面。其缺點是：精度受到空間劃分分辨率的影響，僅適用于小場景，單個物體，遮擋較少的場景，難以處理高精度、大規模的場景。

3.2 基于深度圖融合的方法

基于深度圖融合的方法［43］是利用SFM 產生的稀疏點云構建種子柵格（patch），通過區域生長對patch 進行優化后生成不同視角下的深度圖，再將不同視角深度進行融合的稠密點云重建方法。該方法主要包括以下3 個步驟：數據預處理、區域生長、深度圖融合，如圖15所示。

Fig.15 Dense reconstruction based on depth map fusion圖15 基于深度圖融合的稠密重建

3.2.1 數據預處理

patch 是以空間中的3D 點中心建立一個很小的柵格，柵格的分辨率與圖像分辨率以及3D 點的位置相關。patch的中心點即為3D 點的位置，patch 的朝向即為3D 點的法向量，patch 上的每個點可以投影到不同視角中計算NCC 用于度量光度一致性。

其預處理包括：通過SFM 重建的稀疏點構建種子patch、構建不同分辨率的圖像以處理圖像分辨率不同情況、全局視角選擇。全局視角選擇是指選擇滿足特定條件的視角用于區域生長。一般要滿足以下條件：①圖像具有相同的內容，外觀——通過共享sift 特征點的個數；②圖像具有足夠大的視差（寬基線）——通過視線的夾角；③圖像具有相似的分辨率——通過計算投影球半徑。

為使得計算加速，在生成某一視角的深度圖時，文獻［43］從以上全局視角中選擇4 個視角，稱為局部視角。局部視角的選擇要求：NCC值大于一定的閾值，并且和已經選擇的視角的極平面要足夠分散（不共面），如圖16所示。

Fig.16 Diagram of view selection圖16 視角選擇示意圖

3.2.2 區域生長

在建立初始種子patch 后，可以利用種子patch 逐步擴張的方法在其鄰域生成新的patch。其主要步驟包括：①對稀疏點云中的種子patch 以置信度（光度一致性）建立優先級隊列；②估計初始稀疏種子點的深度；③對每個種子點進行非線性深度優化；④每次優化完后，如果圖像上的對應像素沒有深度值或當前像素的置信度值高于鄰域像素一定范圍，則將其像素添加到隊列中。

其中，非線性優化的數學模型為：

其中，i、j是patch 中的采樣點，k為視角個數。IR(i，j)表示點在參考圖像上的光度，ck是顏色尺度，由于IK(i，j)表達式中含有深度信息，因而可以利用梯度下降等優化方法對E 進行優化，進而求解像素點深度。優化求解結果如圖17所示［43］。

Fig.17 Reconstruction effect based on depth map fusion圖17 基于深度圖融合的重建效果

3.2.3 深度圖融合

深度融合是指在三維空間合并位置相近的點。融合過程同樣要對點施加一定的約束，通常包括一致性約束、可視性約束。一致性約束是指相同的點在不同視角下的深度應當具有一致性；可視性約束是指在圖像上可視的三維點在融合后的深度圖上不能被遮擋。如圖18 所示，A 和C 將被認為是錯誤的點，予以剔除，B 將被選擇為正確的點。

Fig.18 Diagram of consistent point selection圖18 一致點選擇示意圖

深度圖融合方法的特點為：鄰域視角選擇提升了深度估計的準確度，對一些有遮擋、障礙物較多的場景能較好地進行處理，適用場景廣泛，只用到光度一致性約束和可視性約束，算法實現簡單。

3.3 基于空間patch擴散的方法

與深度圖中在二維圖像上進行patch 擴張不同的是，基于空間patch 擴散的方法［44］采用的是在空間中對三維patch（大小為5x5 共25 個三維點）進行擴散，最終使patch覆蓋物體表面。

重建流程如圖19所示。

Fig.19 Reconstruction process based on patch diffusion圖19 基于空間patch擴散重建流程

初始種子點生成采用SIFT 等特征點；擴張過程對已重建三維點的鄰域進行匹配，濾波過程采用兩種約束去除噪聲點：光度一致性約束和可視性約束。

初始3D patch 的生成步驟：①在圖像上均勻計算SIFT/Harris 特征；②沿極線進行搜索找到匹配特征點；③對匹配對，通過三角化建立patch，法向量指向參考圖像，通過光度一致性約束對可視圖像進行篩選；④對patch 位置和法向量進行優化。

空間patch 擴散步驟為：①將三維patch 投影到圖像上；②如果相鄰cell 沒有patch 且深度連續，則進行patch 擴散，新建patch 的法向量初始值等于領域patch，新建patch的位置設置為當前cell 的視線和鄰域patch 所在平面的交點；③計算初始patch 的可視圖像，并進行優化。

擴張過程中并沒有考慮到可視性的情況，濾波過程的主要目的是利用各種約束對這些不合理的點進行篩選。可用的約束包括可視性約束，可視圖像個數小于一定閾值，圖像鄰域中的cell同時也是空間鄰域，其比例小于一定閾值（0.25）。對于每一個patch，統計它所在及相鄰image cell 中的所有patch。如果3D 空間中與相鄰的patch 比例小于一定值（0.25），則認為是孤立點，將被移除。

PMVS 和基于深度圖的重建方法類似，基本原理也是基于光度一致性，并且采用區域生長的方式進行稠密重建，同樣適用于雜亂有遮擋的大規模場景。算法優點是：適用性強，適用于各種形狀的物體，能夠直接重建出完整的稠密點云，不需要進行單獨的深度圖融合；重建過程中同時考慮了可視性約束，能夠得到更加準確的重建結果。其缺點是不適用于非朗伯面（Non-Lambertian），并容易產生空洞，計算量非常大。

近年來，也有針對以上算法不足的改進方法，如：文獻［35］利用并行計算提升點云重建速度；文獻［45］針對在小范圍場景三維重建過程中存在離群點的現象，將PMVS 算法與統計分析法相融合提出一種改進的點云濾波算法，提高了重建精度。

當前，稠密點云重建面臨的主要問題與挑戰是：①弱紋理和朗伯面——不能滿足光度一致性約束；通常需要先驗形狀約束進行重建；②細長結構重建——很難得到像素級別精度的深度值；通常的解決方法需要進行語義級別的重建；③動態場景重建——圖像序列的重建方式可以捕捉動態場景，但是基本只適合于實驗室環境；動態物體的重建需要結合形狀先驗、高級語義信息等。

4 三維表面重建

在得到稠密點云后，利用各種算法可以從點云重建物體表面。三維物體的表面表達方式包括：邊界表示法、空間劃分法、構造體素法。邊界表示法用四邊形、三角形、參數曲面或者非參數的曲面表示物體形狀，具有穩定性強、靈活性強、有助于表示表面細節。空間劃分法將物體內部劃分成細小連續實體以描述物體形狀。構造體素法通過集合運算（并、交、差）將簡單形狀組合成復雜形狀，只能用來表示較為簡單的實體，無法用于形狀復雜或表面精細的物體。

計算機視覺中最常用的三維模型表示方式為三角網格（Triangular Mesh）。其基本結構包括：頂點（Vertex）、面片（Facet）、邊（Edge）。點、線和面上的各種屬性包括：顏色、法向量、紋理坐標（Texture Coordinate）。

常用的三維表面重建方法有：基于二元分割的表面重建方法、基于符號距離場的表面重建方法。

4.1 基于二元分割的表面重建方法

其主要思想是將空間劃分成一組小單元（體素/四面體），采用二元分割思想將這些小單元劃分成內部和外部兩類，介于內部和外部之間的面即為物體表面，如圖20所示［43］。

Fig.20 Surface reconstruction based on binary segmentation圖20 基于二元分割的表面重建

在稠密點云的基礎上進行表面重建主要經歷的步驟如下：

4.1.1 德勞內三角剖分

德勞內三角剖分（Delaunay Triangulation）是指將點云中的相鄰點用一種特定的方式連線，使得任意相鄰的4 個點構成空間中的一個四面體。德勞內三角剖分具有以下特性：

（1）空球特性。任意4 個點的外接球內都不含有其他點。

（2）最大化最小角。德勞內三角剖分最大化四面體上三角形的最小角。

（3）對偶性。德勞內三角剖分和維諾圖（Voronoi Diagram）是互為對偶的關系，如圖21所示［43］。

Fig.21 Duality of Delaunay triangulation and Voronoi diagram圖21 德勞內三角剖分與維諾圖的對偶關系

文獻［46］給出了一種增量式進行德勞內三角剖分的方法。

三維表面的生成問題將轉化為四面體的二元標記問題，即將一剖分四面體標記為“外部”，另一剖分四面體標記為“內部”，內外交界處的空間曲面即為物體表面。

4.1.2 用圖分割方法提取物體表面

圖分割（Graph Cut）是指對一特殊的有向加權圖的某些邊進行切分，將圖分割為沒有交集的兩個部分，使得邊的權重損失最小［46-48］。這一圖的特殊之處在于包含兩個特殊的節點：源節點和匯節點，其中源節點只有出邊，而匯節點只有入邊，如圖22 所示，邊的粗細表示權重大小，切分目標是選取權重最小的邊切分。

Fig.22 Diagram of graph cut圖22 圖分割示意圖

在對稠密點云進行德勞內三角剖分后，根據其與維諾圖的對偶關系生成維諾圖：即使四面體的中心作為維諾圖各小區域的中心，四面體的三角形表面作為維諾圖的邊。在一定約束條件下，對維諾圖的邊賦予一定權重并將其看作有向加權圖，于是便可進行圖分割提取物體表面，如圖23所示［43］。

Fig.23 Surface reconstruction based on graph cut圖23 基于圖分割的表面重建

對維諾圖的邊賦予權重的方法包括以下幾種：可見性約束、光度一致性約束、平滑性約束［46］、剪影約束［43］、表面質量約束［47］。

重建結果如圖23所示［43］。

4.1.3 網格細節優化

用二元分割方法提取的物體表面，由于點云存在噪聲或者空洞，會導致原始網格存在噪聲或者空洞，無法捕捉場景細節，因而需要對網格細節進行優化（Mesh Refinement）。

文獻［44］提出一種基于光度一致性的網格細節優化方法，即利用多視角圖像再次計算相同點在不同視角下的重投影誤差，通過同時或者逐個移動網格頂點的位置，使重投影誤差最小，也使光度一致性最好。度量光度一致性的能量函數為：

其中，h(Ii，)(x)表示圖像i和j之間在像素x處與光度一致性呈單調遞減的某一函數；表示將圖像i通過曲面S重投影到圖像i上；表示重投影的有效區域；Eerror(S)的含義即為所有兩兩視角對下相同區域的重投影誤差之和；Eerror(S)越小，說明生成的網格越準確。

要優化網格頂點位置，即要求以上能量函數關于網格頂點的梯度，應利用梯度下降法進行優化。但由于網格頂點的位置變化實際上會影響其周圍三角面片上點的位置，于是先計算關于某頂點周圍三角面片點的梯度，再采用插值方法計算這一頂點的梯度，插值方法采用重心坐標法。Eerror（S）對［Xi］的梯度為其對網格上Xi一周所有點梯度的加權和具體如圖24所示。

Fig.24 Grid vertex position optimization圖24 網格頂點位置優化

其中，v(x)表示頂點Xi一周三角網格上的任意一點，φi(x)表示v(x)相對于Xi的重點坐標，?E是Eerror(S)相對于v(x)的梯度。

基于二元分割的表面重建方法具有以下特點：①算法流程簡單，容易實現；②適用于大規模場景以及小物體重建，能夠處理復雜表面；③重建效果依賴德勞內三角剖分的質量，容易受到噪聲和外點影響；④通常需要結合網格細節優化獲取更加精細的表面細節。

4.2 基于符號距離場的表面重建方法

基于符號距離場的表面重建流程如圖25所示。

Fig.25 Surface reconstruction based on symbol distance field圖25 基于符號距離場的表面重建流程

Fig.26 Uniform（left）and non-uniform（right）spatial division圖26 均勻（左）與非均勻（右）的空間劃分

4.2.1 空間劃分

空間劃分主要是用于對空間中的隱函數進行離散化，以便于對劃分區域的每個頂點計算一個符號距離值。空間劃分方式可以采用均勻和非均勻兩種。均勻劃分是將空間劃分為大小相等的體素網格（Grid），非均勻方式是用八叉樹（Octree）結構將空間進行劃分，即在點密集處空間劃分較細，樹的頂點較小；點稀疏處空間分辨率較低，樹的頂點較大，點的分辨率可以根據三維點到相鄰點的平均距離確定。

4.2.2 符號距離場構建

文獻［49］提出一種局部的符號距離場構建方法，用符號距離函數（Signed Distance Function，SDF）計算空間中各處的符號距離值，然后提取等值面作為物體表面。

符號距離函數是某度量空間中點X 到某一邊界的距離關于X 位置的函數，當X 在邊界內時，SDF 為正；當X 在邊界外時，SDF 為負。如圖27所示，0所在位置即可認為是物體表面。

Fig.27 Diagram of symbol distance function圖27 符號距離函數示意圖

文獻［49］所用的符號距離函數為：

即在點x附近鄰域取i個點pi，在pi處建立局部坐標系（x軸與pi點的法向量ni方向相同，y軸與z軸分別與x軸垂直），再利用旋轉矩陣Ri將x旋轉到局部坐標系下，得到xi=Ri·(x-pi)，對以下基函數f(xi)用ciwi(xi)加權計算即可得到點x處的符號距離值，其中，wi(xi)權重函數，ci為點pi的置信度（NCC 值）。

局部符號距離場結果如圖28所示［49］。

Fig.28 Boundaries computed by symbolic distance function圖28 符號距離函數計算邊界

文獻［50］提取了一種全局的符號距離場構建方法，將有向點看作是對隱函數梯度的采樣，如圖29 所示。全局的方法優點是魯棒，能夠很好地處理噪聲和空洞，缺點是計算量大，需要求解大規模的方程。

Fig.29 Poisson surface reconstruction圖29 泊松表面重建

4.2.3 移動立方體算法生成表面

移動立方體算法（Marching Cube）最早由文獻［51］提出，在得到每個體素8 個頂點的符號距離值后，通過插值方法生成三角形面片，每個頂點有兩種狀態，總共有256種，但由于反轉狀態不變，因此可以減少一半，為128 種，再根據旋轉不變形，又可以減少到14 種情況。可以認為這14 種類似于基，經過旋轉、反轉可以得到256 種狀態所對應的結果（其中3種情況），如圖30所示。

Fig.30 Surfaces generated by marching cube method圖30 移動立方體法生成表面

文獻［52］提出一種專門針對八叉樹的方法，旨在解決相鄰節點分辨率不一致時可能出現的裂縫問題。通過引入一組從八叉樹的拓撲推導而來的二叉邊樹，可以在不約束八叉樹的拓撲結構或修改頂點值的情況下，直接提取水密的多邊形網格。

文獻［53］提出一種將平面投影與區域生長相結合的表面重建方法，降低了表面重建復雜度。文獻［54］用統計濾波器對點云簡化去噪，并建立點云間拓撲結構，緩解了重建表面粗糙、孔洞等問題。

5 紋理圖像創建

為了生成更加逼真的場景和物體，可以對生成的三維模型貼上紋理圖像。文獻［56］在文獻［55］的基礎上提出創建紋理圖像管線，主要分為如下步驟：視角選擇、全局顏色調整、泊松圖像編輯。

5.1 視角選擇

為每個三角面片選擇唯一的視角，用于獲取紋理信息，視角的選擇應當考慮到以下因素：圖像尺度、圖像細節豐富程度、圖像可視性、鄰域平滑性。具體方法為：將視角選擇建模成多標簽分配問題，為每一個三角形分配一個標簽（視角），在網格上建立馬爾可夫隨機場，每個面片代表一個頂點，通過優化能量函數得到最優視角選擇。能量函數分為數據項Edata和平滑項Esmooth，即：

其中，Fi表示第i個三角面片，li表示分配給它的視角標簽，數據項Edata通過三角面片對應的二維圖像上三角形內的平均像素梯度進行度量，旨在選擇圖像分辨率較大的視角；平滑項約束希望相鄰三角形有相同的視角以保持鄰域的平滑性，從而減少縫隙，為后續處理提供便利。選擇好對應的視角后，即可將網格投影到對應視角的可視圖像上，截取對應圖像作為紋理圖像，經過坐標歸一化處理后，即可作為三角面片的紋理坐標，后期通過紋理坐標可找到紋理圖像上對應的貼片。

5.2 全局顏色調整

由于不同視角間存在相機曝光或者光照差異導致不同的紋理圖像邊界處存在明顯縫隙，因而要為邊界處的每個像素添加一個顏色調整量使得縫隙處的顏色差異盡量小，紋理圖像內部相鄰像素的調整量盡量相似。

具體方法：首先為三角面片的每個頂點對應像素添加一個調整量，再通過插值方式得到三角面片內部對應的每個像素的調整量，然后將縫隙處的頂點拆分成n個（n為標簽個數），對其施加n個調整量，最后通過優化以下能量函數得到調整后的結果：

其中，Es(g)的約束目標是使不同視角下同一頂點調整后的顏色盡量接近，而Ei(g)的約束目標是使相同視角下不同頂點的調整量盡量接近，λ為可調整系數。

5.3 泊松圖像編輯

對于縫隙比較嚴重的區域，全局顏色調整并不能保證完全去除縫隙，還需要作進一步處理。泊松圖像編輯算法原理為對前景圖像和背景圖像進行混合，使得融合后的圖像滿足：邊界上的像素值與背景圖像相同，前景區域內的梯度與引導梯度場相同［57］。

這種方法最初用于圖像融合。如圖31 所示［57］，要將源圖像g插入到目標圖像S的區域Ω中，該區域具有邊界?Ω，目標圖像的像素值與像素坐標的關系可表示為一個函數f*，而合成后的圖像在區域Ω中的像素值與像素坐標關系用f描述，v代表一個引導梯度場，實際上是指源圖像g的梯度。其中，g、S、Ω、?Ω、f*、v都是已知量，而f是待求函數（融合后在區域Ω 中的圖像）。

Fig.31 Poisson image fusion圖31 泊松圖像融合示意圖

優化的函數為：

其中，?f表示圖像融合后f的梯度，優化目標是使f邊界上的像素值與背景圖像相同，f的梯度與源圖像g的梯度v相同。求解以上優化函數可以轉化為求解以下泊松方程：

其中，Δf表示f的拉普拉斯濾波結果，在某一像素處可以用相鄰像素線性表示，而divv==Δg，于是可將泊松方程轉化成Δf=Δg，用線性方程快速求解。泊松圖像融合效果如圖32所示［57］。

Fig.32 Effect of Poisson image fusion圖32 泊松圖像融合效果

進行全局顏色調整后，只需在縫隙處進行泊松編輯即可得到邊界平滑的紋理圖像。此外，可以通過控制內部區域的范圍（如由當前邊界往前景區域擴充3 個像素作為前景進行計算），從而減少計算復雜度。泊松圖像編輯還應用于圖像拼接，如文獻［58］提出一種對無人機影像的拼接技術，取得了良好效果。

6 結語

與深度學習的三維重建方法相比，基于圖像的三維重建算法流程清晰，具有明確的中間結果，是典型的“白盒”模型，不需要預訓練神經網絡，對各種場景和物體適應性強，可以全流程優化改善，增強實現效果。但也存在計算速度慢、過分依賴于圖像特征點，且對弱紋理、朗伯面、纖細物體重建效果不佳等缺點。

鑒于以上不足，以下工作可能取得重要進展：①圖像數據與視覺里程計、激光雷達等多傳感器信息融合，不再依賴于純圖像，增加信息來源，既提升運算速度，又提高計算精度；②與深度學習方法相結合，實現在沒有圖像特征點或匹配對很少的情況下，得到相機的相對位置關系，并對弱紋理區域、朗伯面重建具有魯棒性；③算法軟件提升與基于GPU 的硬件加速相結合，提升算法運行效率，逐步實現大場景下的實時三維重建。