多視圖的三維景物中平表面重建

李靜，楊宜民，蔡述庭

(1.洛陽師范學(xué)院 信息技術(shù)學(xué)院，河南 洛陽 471022； 2. 廣東工業(yè)大學(xué) 自動化學(xué)院，廣東 廣州 510090)

重建位于同一張平表面上的三維空間點(diǎn)，這對于室內(nèi)場景以及一些建筑場景是很常見的情形。傳統(tǒng)的方法是先重建出三維空間點(diǎn)[1-2]，再用這些三維空間點(diǎn)去擬合場景中的平表面[3-4]。此方法存在2個缺陷：1)重建三維空間點(diǎn)時，沒有利用這些三維空間點(diǎn)位于一個平表面上這一先驗知識；2)在擬合景物平表面時，需測量三維空間點(diǎn)之間距離的大小，在射影空間中，三維空間點(diǎn)之間的距離沒有物理意義。針對這些缺陷，文獻(xiàn)[5]將三維空間點(diǎn)約束到一個已知平表面上，并將問題轉(zhuǎn)化為一個8次多項式進(jìn)行求解。文獻(xiàn)[6]考慮到三維點(diǎn)位于景物平表面的約束條件，針對場景平表面已知與景物平表面未知2種情形，給出了恢復(fù)三維景物平表面的方法，但不能保證獲得全局最優(yōu)。文獻(xiàn)[7]利用場景平表面與單應(yīng)矩陣之間的關(guān)系，通過求解單應(yīng)矩陣來求解場景平面，由于單應(yīng)矩陣的轉(zhuǎn)移誤差函數(shù)是擬凸函數(shù)，可以利用最小化L范數(shù)方法估計單應(yīng)矩陣，繼而恢復(fù)景物平表面，利用L方法能夠得到全局最優(yōu)值，缺點(diǎn)是計算效率較低，且對局外點(diǎn)敏感。另外，文獻(xiàn)[5-7]都是針對2幅圖像的三維場景平表面重建，且不能推廣到多幅圖像。因為考慮三維點(diǎn)位于一張已知的景物平表面上，則對應(yīng)點(diǎn)集的反向投影線不僅要相交，而且要相交于一個平表面上，以此為約束條件，給出最小化反投影誤差的重建法和最小化轉(zhuǎn)移誤差的重建法，并采用GA算法進(jìn)行了最小化求解。采用最小化反投影誤差的平表面重建法記為“RPE-GA”(reverse projection error-GA)法，采用最小化轉(zhuǎn)移誤差的平表面重建法記為“TE-GA”(transfer error-GA)法，實(shí)驗結(jié)果證明了所提方法的優(yōu)越性。

1 基于兩視圖的三維景物中平表面重建

如圖1所示，一張景物平表面上的點(diǎn)的圖像點(diǎn)與在第2幅視圖上的對應(yīng)圖像點(diǎn)由一個單應(yīng)矩陣相關(guān)聯(lián)，事實(shí)上，這一景物平表面誘導(dǎo)了2幅視圖之間的一個單應(yīng)矩陣。來自空間不同景物平表面的圖像對應(yīng)點(diǎn)之間滿足不同的平面單應(yīng)約束。

圖1 由景物平表面誘導(dǎo)的單應(yīng)矩陣H(π)Fig.1 The homography H(π) induced by a scene plane

是空間平表面上的點(diǎn)在2幅圖像上的對應(yīng)點(diǎn)對，π為未知景物平表面，空間景物平表面參數(shù)化為

π=[nTd]T

式中：n∈R3(法向量)，0≠d∈R。則存在矩陣H(π)使得

x′i?H(π)xi

式中:?表示相差一個比例因子的情況下相等，H(π)是一個3×3的矩陣，表示從圖像1到圖像2通過任意平面π的轉(zhuǎn)移映射。同理，從圖像2到圖像1的轉(zhuǎn)移映射表示：

xi?H-1(π)x′i

若兩視圖投影矩陣為

P1=K[I0]，P2=K′[Rt]

則單應(yīng)矩陣H(π)表示為

H(π)?K′(R-t(n/d)T)K-1

證明為計算H(π)，把第1幅視圖上的點(diǎn)反向投影且確定該射線和平面π=[nTd]T的交點(diǎn)X，然后再把這個3-D點(diǎn)X投影到第2幅視圖上。對第1幅視圖有x?P1X=K[I0]X，因此該射線上的任何點(diǎn)x=[xTK-Tρ]T都可以投影到x，其中ρ是確定該射線上某個點(diǎn)的參數(shù)，因3-D點(diǎn)X在平面π上，所以滿足πTx=0，則ρ=-(nTK-1x)/d，從而，

3-D點(diǎn)X投影到第2幅圖像上可得到

由此，可得H(π)?K′(R-(tnT)/d)K-1x。

若不存在噪聲干擾和誤匹配，一般通過最小化幾何誤差來找到一個單應(yīng)矩陣H(π)，即可求得空間平面參數(shù)π，若

則最小化幾何誤差函數(shù)表示為

這種方法的目標(biāo)函數(shù)為L2范數(shù)，所以稱之為L2方法，缺陷是容易陷入局部最優(yōu)。將求解單應(yīng)矩陣的問題轉(zhuǎn)換為最小化L范數(shù)進(jìn)行求解，即最小化最大轉(zhuǎn)移誤差來求解，表示為

min maxei(π)

式中：ei(π)是轉(zhuǎn)移誤差函數(shù)，可表示為

(1)

2 基于多視圖的三維景物中平表面重建

三維景物中平表面重建意味著重建后的空間點(diǎn)應(yīng)位于同一個平面內(nèi)。換言之，三維景物中平表面重建就是以重建的空間點(diǎn)在同一個平面內(nèi)為約束條件的尋優(yōu)問題。采用基于最小化反投影誤差的平表面重建法記為“RPE-GA”，采用基于最小化轉(zhuǎn)移誤差的平表面重建法記為“TE-GA”。

2.1 基于最小化反投影誤差的平表面重建模型

(2)

由于在圖像點(diǎn)檢測與匹配的過程中，不可避免地存在誤差，所以這些反投影直線并不相交于一點(diǎn)，如圖2所示。一般的求解方法是通過最小化點(diǎn)到所有反投影直線的距離平方和來得到最優(yōu)解Xj，但這種方法只適用于歐氏空間或者準(zhǔn)歐氏空間，如果未得到攝像機(jī)的標(biāo)定參數(shù)，即僅僅在射影空間中進(jìn)行三維重建，此時歐氏距離并沒有實(shí)際的物理含義，因此不能保證此方法的最優(yōu)性。因此很多文獻(xiàn)[9-12]都是將誤差轉(zhuǎn)移到圖像平面上(即歐氏空間中)進(jìn)行最小化重投影誤差來求解的。本文是針對景物中的平表面進(jìn)行重建，由于誤差的存在，反向投影線雖然也不交于一點(diǎn)，但反向投影線都應(yīng)與一個空間平面相交，如圖3所示。因為反投影線都與一個場景平面π相交，且所求空間點(diǎn)也位于場景平面π上，因此平表面重建可以通過最小化平表面上的反投影誤差進(jìn)行求解，表示為

i=1,2,…,M;j=1,2,…,N

(3)

圖2 反投影線不交于一點(diǎn)Fig.2 Reverse projection lines don’t intersect in a point

圖3 反向投影線與一個空間平面相交Fig.3 Reverse projection lines intersect in a scene plane

2.2 基于最小化轉(zhuǎn)移誤差的平表面重建模型

因此，平面π到第1個圖像平面的單應(yīng)矩陣為H1=K1(R1-t1(n/d)T)。同理，在第2個攝像機(jī)下的投影表示為

i=1,2,…,M;j=1,2,…,N

(4)

由于可以很明顯地看出，式(4)較式(3)少了未知參數(shù)αij，求解相對簡單一些。式(4)是利用二維圖像平面與二維空間平面之間的轉(zhuǎn)移矩陣(即單應(yīng)矩陣Hi(π))將誤差轉(zhuǎn)移到空間平表面上進(jìn)行最小化的。而式(3)是通過限制反投影線與空間平面相交且交于一點(diǎn)，從而將誤差轉(zhuǎn)移到空間平面上最小化求解。總而言之，無論最小化反投影誤差的平表面重建法，還是最小化轉(zhuǎn)移誤差的平表面重建法，都是將誤差轉(zhuǎn)移到空間平表面上進(jìn)行最小化求解的，這2種方法的基本原理一致。

2.3 基于多視圖的三維景物中平表面重建

式(3)和式(4)都是有約束的優(yōu)化問題，對于有約束的最小化問題，可以引入罰函數(shù)構(gòu)造增廣代價函數(shù)，將給定的約束優(yōu)化問題轉(zhuǎn)化為一系列無約束優(yōu)化問題，然后通過求解這一系列的無約束優(yōu)化問題而獲得原約束優(yōu)化問題的解。懲罰法有多種類型，常用的有：外懲罰法(又稱罰函數(shù)法)和內(nèi)懲罰法(又稱障礙函數(shù)法)。采用罰函數(shù)法，如果約束優(yōu)化問題為

minf(x)

s.t.gi(x)=0,i=1,2,…,M

s.t.hij(x)=0,i=1,2,…,M;j=1,2,…,N

(5)

記可行域為Ω={x∈Rn|gi(x)=0,hij(x)=0,i=1,2,..,M;j=1,2,…,N}，首先選擇一個充分大的數(shù)β，構(gòu)造一個懲罰函數(shù)：

式中:β為懲罰因子，表示對不可行點(diǎn)的懲罰。然后，將原問題轉(zhuǎn)換為求解下列無約束優(yōu)化問題：

minb(x)

(6)

稱式(6)為式(5)的增廣代價函數(shù)。由于β是一個非常大的數(shù)，所以函數(shù)b(x)的解也一定是原問題的解。

先將式(3)和式(4)轉(zhuǎn)換為無約束問題，再采用遺傳算法(GA)進(jìn)行求解，遺傳算法是一個非常成熟的計算工具，這里不再詳細(xì)介紹。本文的創(chuàng)新之處是建立了2種平表面重建模型，而不是針對遺傳算法的改進(jìn)。另外，由法向量n中的3個參數(shù)足以確定一個平面，所以在實(shí)際計算過程中令d=1。

基于多視圖的三維景物中平表面重建步驟為：

1)輸入N個空間點(diǎn)對應(yīng)的圖像點(diǎn)數(shù)據(jù)，攝像機(jī)的投影矩陣Pi,i=1,2,…,M；

2)將式(4)和式(3)轉(zhuǎn)換為如式(6)所示的無約束優(yōu)化函數(shù)，設(shè)置懲罰因子β=500，即構(gòu)造關(guān)于式(4)和式(3)的適應(yīng)度函數(shù)；

3)初始化種群，每一個染色體都是N個空間點(diǎn)和場景平面法線向量n，以及參數(shù)αij(見式(3))的組合，采用實(shí)數(shù)編碼，設(shè)定種群個數(shù)，即染色體個數(shù)Nind=50，并設(shè)置邊界條件；

4)選擇操作，計算種群中每個染色體的適應(yīng)度值，采用輪盤賭的方法進(jìn)行選擇復(fù)制，即各個染色體按照適應(yīng)度值所占總適應(yīng)度值的比例組成面積為1的一個圓盤，指針轉(zhuǎn)動停止之后，指向的染色體將被復(fù)制到下一代，適應(yīng)度好的染色體被選擇的概率大，每個染色體被選中的概率為

式中：f(xi)代表適應(yīng)度值，Nind代表染色體個數(shù)；

5)交叉操作，以交叉概率pc=0.6進(jìn)行一點(diǎn)交叉，把2個父代個體的部分基因進(jìn)行交換而生成新的個體，直到所有個體都被選擇過；

6)變異操作，選變異概率pm=0.02，即全部染色體位串上的每位基因按變異概率pm進(jìn)行隨機(jī)改變；

7)保存當(dāng)前適應(yīng)度函數(shù)的最小值，以及對應(yīng)的染色體位置；

8)直至達(dá)到最大迭代次數(shù)，最大迭代次數(shù)為100，若滿足，則停止迭代，否則轉(zhuǎn)4)；

3 實(shí)驗與結(jié)果分析

采用2.6 GHz Pentium Dual-Core CPU，2 GB內(nèi)存的臺式電腦，并在MATLAB 7.0環(huán)境下進(jìn)行了2組實(shí)驗，分別是仿真數(shù)據(jù)實(shí)驗和真實(shí)圖像實(shí)驗。并對式(3)和式(4)利用GA算法進(jìn)行了優(yōu)化計算，采用式(3)記為“RPE-GA”法，采用式(4)記為“TE-GA”法。為了進(jìn)行對比，利用傳統(tǒng)法進(jìn)行實(shí)驗，首先重建出三維空間點(diǎn)，然后擬合一個平表面，在傳統(tǒng)法中利用L方法[11]進(jìn)行多視圖的重建，L方法可以獲得全局最優(yōu)解，因此將這種方法記為“L”法。實(shí)驗中采用的性能參數(shù)為：重投影均方根誤差為i=1,2,…,M，j=1,2,…,N，xij代表檢測點(diǎn)坐標(biāo)，代表反投影點(diǎn)坐標(biāo)，單位為像素；空間位置均方根誤差為

3.1 仿真實(shí)驗與結(jié)果分析

仿真實(shí)驗中需要知道產(chǎn)生檢測數(shù)據(jù)的真正平表面參數(shù)，設(shè)置如下：將平表面定為Xw-Yw平表面，即Zw=0平表面，隨機(jī)生成范圍在[-1,1]內(nèi)的100個點(diǎn)，即N=100，要保證這些點(diǎn)在每個攝像機(jī)下都是可視的。假設(shè)用6個攝像機(jī)進(jìn)行拍攝，實(shí)驗中每臺攝像機(jī)的內(nèi)參數(shù)設(shè)為

設(shè)定每個攝像機(jī)之間大概相差10°，圖像分辨率為640×480。仿真實(shí)驗環(huán)境如圖4所示，圖4中用圓圈來表示攝像機(jī)的位置，點(diǎn)云為隨機(jī)生成的Zw=0平表面上的點(diǎn)。

圖4 仿真實(shí)驗環(huán)境示意圖Fig.4 The schematic diagram of simulation environment

在每個圖像點(diǎn)上增加零均值、標(biāo)準(zhǔn)差為σ的高斯噪聲，噪聲大小分別為0.2、0.4、0.6、0.8、1.0個像素，對每種噪聲分別進(jìn)行100次實(shí)驗。由于對帶有噪聲的數(shù)據(jù)都進(jìn)行了100次實(shí)驗，每次實(shí)驗都會產(chǎn)生一個重建結(jié)果，數(shù)據(jù)過多，不再列出，只給出3種平表面重建方法的性能比較，如表1所示。從表1可以看出，RPE-GA法和TE-GA法的精度基本一致，可以明顯看出空間位置均方根誤差比較小，這是因為是通過最小化的空間點(diǎn)位置誤差進(jìn)行優(yōu)化求解的，所以空間位置均方根誤差較小；而L方法的重投影均方根誤差較小，這是因為在L法中是通過最小化圖像平面上的重投影誤差進(jìn)行優(yōu)化求解的，所以重投影均方根誤差較小。由于在L法中沒有考慮到景物平表面的約束，所以空間位置均方根誤差較大，也就與真實(shí)的景物平表面相差較遠(yuǎn)。

表1 3種平表面重建方法的性能比較

3.2 真實(shí)圖像實(shí)驗與結(jié)果分析

真實(shí)圖像實(shí)驗采用Zhang從不同角度拍攝平面模版的5幅圖像(http://research.microsoft.com/en-us/um/people/zhang/calib/)，如圖5所示，圖像分辨率大小為640×480，平面模版上共有256個角點(diǎn)，即N=256。

圖5 平面模版的5幅圖像Fig.5 Five images of a planer pattern

真實(shí)圖像的重建結(jié)果如圖6所示。從圖6可以看出本文方法的重建結(jié)果基本上位于一個平表面上，而傳統(tǒng)的L法重建結(jié)果非常散亂，并不位于一個平表面上。真實(shí)圖像實(shí)驗中3種重建平表面方法的性能比較如表2所示。同樣的，RPE-GA法和TE-GA法的空間位置均方根誤差比較小，與真實(shí)的景物平表面相接近；而L法的重投影均方根誤差較小，但空間位置均方根誤差卻比較大，也就是與真實(shí)的景物平表面相差較遠(yuǎn)。從原理上來講，RPE-GA和TE-GA 2種優(yōu)化方法是一樣的，都是將空間點(diǎn)限制在一個平表面上，然后最小化空間點(diǎn)之間的距離誤差，所以它們的精度基本一致。只是由于RPE-GA法中多了參數(shù)αij，所以RPE-GA法的計算會比TE-GA法復(fù)雜。本文方法的優(yōu)越性：1)利用三維空間點(diǎn)都位于同一個景物平表面這一先驗知識，計算出的平表面精度較高；2)可以推廣到多幅視圖。

需要指出的是RPE-GA法和TE-GA法對局外點(diǎn)都比較敏感，事實(shí)上，目前已有的基于兩視圖景物平表面重建方法，如Olsson[7]提出的重建方法對局外點(diǎn)也比較敏感。如果用在場景比較復(fù)雜的真實(shí)圖像實(shí)驗中，會存在2類局外點(diǎn)，一類是在匹配錯誤的點(diǎn)，一類是不位于一個平表面上的點(diǎn)。首先需要采用一些剔除局外點(diǎn)的方法對局外點(diǎn)進(jìn)行剔除，Olsson采用L1方法對局外點(diǎn)進(jìn)行了剔除，但是從實(shí)驗結(jié)果中可以很明顯的看出，即使對局外點(diǎn)進(jìn)行了剔除，仍然有很多圖像點(diǎn)不位于同一個平表面上，用這些不位于同一個平表面上的圖像點(diǎn)去重建一個平表面勢必是不準(zhǔn)確的。

(a) L的重建結(jié)果

(b) RPE-GA的重建結(jié)果

(c) TE-GA的重建結(jié)果 圖6 平面模版的重建結(jié)果Fig.6 Reconstruction of the planar pattern

Table2Theperformanceofthreeplanereconstructionmethods

方法名稱重投影均方根誤差空間位置均方根誤差L挰0.063 50.058 3RPE-GA0.081 20.001 3TE-GA0.081 50.001 2

5 結(jié)束語

在傳統(tǒng)的三維景物中平表面重建法中，都是先重建出三維點(diǎn)，再用這些三維點(diǎn)去擬合景物平表面，由于沒有利用三維空間點(diǎn)位于一個平表面的先驗知識，所以重建出的三維點(diǎn)并不位于一個平表面上。針對這個問題，本文提出了2種三維景物中平表面重建模型，2種模型都是考慮到重建的空間點(diǎn)應(yīng)位于同一個平面內(nèi)，以此為約束條件分別在空間平面上進(jìn)行最小化反投影誤差和最小化轉(zhuǎn)移誤差，再采用GA算法對2種模型進(jìn)行了優(yōu)化求解，從而獲得平表面重建結(jié)果。實(shí)際上，基于最小化反投影誤差的平表面重建法和基于最小化轉(zhuǎn)移誤差法的平表面重建法基本原理相同，只是計算復(fù)雜度不同。實(shí)驗表明了本文方法的優(yōu)越性。另外如何剔除掉局外點(diǎn)，或者說如何從復(fù)雜場景中提取出位于一個平表面上的圖像點(diǎn)，將RPE-GA法和TE-GA法用于復(fù)雜場景的真實(shí)圖像實(shí)驗上，將是未來需要研究的問題。據(jù)知，目前在國內(nèi)外考慮三維點(diǎn)位于空間景物平面的約束條件，針對多視圖的三維景物中平表面重建的研究還處于空白。

參考文獻(xiàn):

[1]LINDSTROM P. Triangulation made easy[C]//Proceedings of IEEE Conference Computer Vision Pattern Recognition. San Fransico, USA, 2010: 1554-1561.

[2]TOSSAVAINEN T. Approximate and SQP two view triangulation[C]//Proceedings of 10th Asian Conference Computer Vision. Queenstown, New Zealand, 2010, 3: 1303-1316.

[3]HARTLEY R, ZISSERMAN A. Multiple view geometry in computer vision[M]. Cambridge: Cambridge University Press, 2003: 54-57.

[4]ZACH C, SORMAN M, KARNER K. High-performance multi-view reconstruction[C]//Proceedings of the Third International Symposium on 3D Data Processing, Visualization, and Transmission. Chapel Hill, North Carolina, USA, 2006: 113-120.

[5]CHUM O, PAJDLA T, STURM P. The geometric error for homographies[J]. Computer Vision and Image Understanding, 2002, 97(1): 86-102.

[6]KANATANI K, NIITSUMA H. Optimal two-view planar scene triangulation[J]. IPSJ Transactions on Computer Vision and Applications, 2011, 3: 67-79.

[7]OLSSON C, ERIKSSON A. Triangulating a plane[C]//Proceedings of Scandinavian Conference on Image Analysis. Berlin Heidelberg, Germany, 2011: 13-23.

[8]KE Q, KANADA T. Quasiconvex optimization for robust geometric reconstruction[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2007, 29(10): 1834-1847.

[9]AGARWAL S, SNAVELY N, SEITZ S M. Fast algorithms for L-infinity problems in multiview geometry[C]//Proceeding of IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Anchorage, Alaska, USA, 2008: 24-26.

[10]HARTLEY R, KAHL F, OLSSON C, et al. Verifying globle minima forL2minimization problems in multiple view geometry[J]. International Journal of Computer Vision, 2013: 288-304.

[11]KAHL F, HARTLEY R. Multiple-view geometry under the L-infinity norm[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2008, 30(9): 1603-1617.

[12]LEE H, SEO Y, LEE S W. Removing outliers by minimizing the sum of infeasibilities[J]. Image and Vision Computing, 2010, 28: 881-889.