李靜,楊宜民,蔡述庭
(1.洛陽師范學(xué)院 信息技術(shù)學(xué)院,河南 洛陽 471022; 2. 廣東工業(yè)大學(xué) 自動(dòng)化學(xué)院,廣東 廣州 510090)
重建位于同一張平表面上的三維空間點(diǎn),這對(duì)于室內(nèi)場景以及一些建筑場景是很常見的情形。傳統(tǒng)的方法是先重建出三維空間點(diǎn)[1-2],再用這些三維空間點(diǎn)去擬合場景中的平表面[3-4]。此方法存在2個(gè)缺陷:1)重建三維空間點(diǎn)時(shí),沒有利用這些三維空間點(diǎn)位于一個(gè)平表面上這一先驗(yàn)知識(shí);2)在擬合景物平表面時(shí),需測(cè)量三維空間點(diǎn)之間距離的大小,在射影空間中,三維空間點(diǎn)之間的距離沒有物理意義。針對(duì)這些缺陷,文獻(xiàn)[5]將三維空間點(diǎn)約束到一個(gè)已知平表面上,并將問題轉(zhuǎn)化為一個(gè)8次多項(xiàng)式進(jìn)行求解。文獻(xiàn)[6]考慮到三維點(diǎn)位于景物平表面的約束條件,針對(duì)場景平表面已知與景物平表面未知2種情形,給出了恢復(fù)三維景物平表面的方法,但不能保證獲得全局最優(yōu)。文獻(xiàn)[7]利用場景平表面與單應(yīng)矩陣之間的關(guān)系,通過求解單應(yīng)矩陣來求解場景平面,由于單應(yīng)矩陣的轉(zhuǎn)移誤差函數(shù)是擬凸函數(shù),可以利用最小化L范數(shù)方法估計(jì)單應(yīng)矩陣,繼而恢復(fù)景物平表面,利用L方法能夠得到全局最優(yōu)值,缺點(diǎn)是計(jì)算效率較低,且對(duì)局外點(diǎn)敏感。另外,文獻(xiàn)[5-7]都是針對(duì)2幅圖像的三維場景平表面重建,且不能推廣到多幅圖像。因?yàn)榭紤]三維點(diǎn)位于一張已知的景物平表面上,則對(duì)應(yīng)點(diǎn)集的反向投影線不僅要相交,而且要相交于一個(gè)平表面上,以此為約束條件,給出最小化反投影誤差的重建法和最小化轉(zhuǎn)移誤差的重建法,并采用GA算法進(jìn)行了最小化求解。采用最小化反投影誤差的平表面重建法記為“RPE-GA”(reverse projection error-GA)法,采用最小化轉(zhuǎn)移誤差的平表面重建法記為“TE-GA”(transfer error-GA)法,實(shí)驗(yàn)結(jié)果證明了所提方法的優(yōu)越性。
如圖1所示,一張景物平表面上的點(diǎn)的圖像點(diǎn)與在第2幅視圖上的對(duì)應(yīng)圖像點(diǎn)由一個(gè)單應(yīng)矩陣相關(guān)聯(lián),事實(shí)上,這一景物平表面誘導(dǎo)了2幅視圖之間的一個(gè)單應(yīng)矩陣。來自空間不同景物平表面的圖像對(duì)應(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幅圖像上的對(duì)應(yīng)點(diǎn)對(duì),π為未知景物平表面,空間景物平表面參數(shù)化為
π=[nTd]T
式中:n∈R3(法向量),0≠d∈R。則存在矩陣H(π)使得
x′i?H(π)xi
式中:?表示相差一個(gè)比例因子的情況下相等,H(π)是一個(gè)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
證明為計(jì)算H(π),把第1幅視圖上的點(diǎn)反向投影且確定該射線和平面π=[nTd]T的交點(diǎn)X,然后再把這個(gè)3-D點(diǎn)X投影到第2幅視圖上。對(duì)第1幅視圖有x?P1X=K[I0]X,因此該射線上的任何點(diǎn)x=[xTK-Tρ]T都可以投影到x,其中ρ是確定該射線上某個(gè)點(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。
若不存在噪聲干擾和誤匹配,一般通過最小化幾何誤差來找到一個(gè)單應(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)
三維景物中平表面重建意味著重建后的空間點(diǎn)應(yīng)位于同一個(gè)平面內(nèi)。換言之,三維景物中平表面重建就是以重建的空間點(diǎn)在同一個(gè)平面內(nèi)為約束條件的尋優(yōu)問題。采用基于最小化反投影誤差的平表面重建法記為“RPE-GA”,采用基于最小化轉(zhuǎn)移誤差的平表面重建法記為“TE-GA”。
(2)
由于在圖像點(diǎn)檢測(cè)與匹配的過程中,不可避免地存在誤差,所以這些反投影直線并不相交于一點(diǎn),如圖2所示。一般的求解方法是通過最小化點(diǎn)到所有反投影直線的距離平方和來得到最優(yōu)解Xj,但這種方法只適用于歐氏空間或者準(zhǔn)歐氏空間,如果未得到攝像機(jī)的標(biāo)定參數(shù),即僅僅在射影空間中進(jìn)行三維重建,此時(shí)歐氏距離并沒有實(shí)際的物理含義,因此不能保證此方法的最優(yōu)性。因此很多文獻(xiàn)[9-12]都是將誤差轉(zhuǎn)移到圖像平面上(即歐氏空間中)進(jìn)行最小化重投影誤差來求解的。本文是針對(duì)景物中的平表面進(jìn)行重建,由于誤差的存在,反向投影線雖然也不交于一點(diǎn),但反向投影線都應(yīng)與一個(gè)空間平面相交,如圖3所示。因?yàn)榉赐队熬€都與一個(gè)場景平面π相交,且所求空間點(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 反向投影線與一個(gè)空間平面相交Fig.3 Reverse projection lines intersect in a scene plane
因此,平面π到第1個(gè)圖像平面的單應(yīng)矩陣為H1=K1(R1-t1(n/d)T)。同理,在第2個(gè)攝像機(jī)下的投影表示為
i=1,2,…,M;j=1,2,…,N
(4)
由于可以很明顯地看出,式(4)較式(3)少了未知參數(shù)αij,求解相對(duì)簡單一些。式(4)是利用二維圖像平面與二維空間平面之間的轉(zhuǎn)移矩陣(即單應(yīng)矩陣Hi(π))將誤差轉(zhuǎn)移到空間平表面上進(jìn)行最小化的。而式(3)是通過限制反投影線與空間平面相交且交于一點(diǎn),從而將誤差轉(zhuǎn)移到空間平面上最小化求解??偠灾瑹o論最小化反投影誤差的平表面重建法,還是最小化轉(zhuǎn)移誤差的平表面重建法,都是將誤差轉(zhuǎn)移到空間平表面上進(jìn)行最小化求解的,這2種方法的基本原理一致。
式(3)和式(4)都是有約束的優(yōu)化問題,對(duì)于有約束的最小化問題,可以引入罰函數(shù)構(gòu)造增廣代價(jià)函數(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)
記可行域?yàn)棣?{x∈Rn|gi(x)=0,hij(x)=0,i=1,2,..,M;j=1,2,…,N},首先選擇一個(gè)充分大的數(shù)β,構(gòu)造一個(gè)懲罰函數(shù):
式中:β為懲罰因子,表示對(duì)不可行點(diǎn)的懲罰。然后,將原問題轉(zhuǎn)換為求解下列無約束優(yōu)化問題:
minb(x)
(6)
稱式(6)為式(5)的增廣代價(jià)函數(shù)。由于β是一個(gè)非常大的數(shù),所以函數(shù)b(x)的解也一定是原問題的解。
先將式(3)和式(4)轉(zhuǎn)換為無約束問題,再采用遺傳算法(GA)進(jìn)行求解,遺傳算法是一個(gè)非常成熟的計(jì)算工具,這里不再詳細(xì)介紹。本文的創(chuàng)新之處是建立了2種平表面重建模型,而不是針對(duì)遺傳算法的改進(jìn)。另外,由法向量n中的3個(gè)參數(shù)足以確定一個(gè)平面,所以在實(shí)際計(jì)算過程中令d=1。
基于多視圖的三維景物中平表面重建步驟為:
1)輸入N個(gè)空間點(diǎn)對(duì)應(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)初始化種群,每一個(gè)染色體都是N個(gè)空間點(diǎn)和場景平面法線向量n,以及參數(shù)αij(見式(3))的組合,采用實(shí)數(shù)編碼,設(shè)定種群個(gè)數(shù),即染色體個(gè)數(shù)Nind=50,并設(shè)置邊界條件;
4)選擇操作,計(jì)算種群中每個(gè)染色體的適應(yīng)度值,采用輪盤賭的方法進(jìn)行選擇復(fù)制,即各個(gè)染色體按照適應(yīng)度值所占總適應(yīng)度值的比例組成面積為1的一個(gè)圓盤,指針轉(zhuǎn)動(dòng)停止之后,指向的染色體將被復(fù)制到下一代,適應(yīng)度好的染色體被選擇的概率大,每個(gè)染色體被選中的概率為
式中:f(xi)代表適應(yīng)度值,Nind代表染色體個(gè)數(shù);
5)交叉操作,以交叉概率pc=0.6進(jìn)行一點(diǎn)交叉,把2個(gè)父代個(gè)體的部分基因進(jìn)行交換而生成新的個(gè)體,直到所有個(gè)體都被選擇過;
6)變異操作,選變異概率pm=0.02,即全部染色體位串上的每位基因按變異概率pm進(jìn)行隨機(jī)改變;
7)保存當(dāng)前適應(yīng)度函數(shù)的最小值,以及對(duì)應(yīng)的染色體位置;
8)直至達(dá)到最大迭代次數(shù),最大迭代次數(shù)為100,若滿足,則停止迭代,否則轉(zhuǎn)4);
采用2.6 GHz Pentium Dual-Core CPU,2 GB內(nèi)存的臺(tái)式電腦,并在MATLAB 7.0環(huán)境下進(jìn)行了2組實(shí)驗(yàn),分別是仿真數(shù)據(jù)實(shí)驗(yàn)和真實(shí)圖像實(shí)驗(yàn)。并對(duì)式(3)和式(4)利用GA算法進(jìn)行了優(yōu)化計(jì)算,采用式(3)記為“RPE-GA”法,采用式(4)記為“TE-GA”法。為了進(jìn)行對(duì)比,利用傳統(tǒng)法進(jìn)行實(shí)驗(yàn),首先重建出三維空間點(diǎn),然后擬合一個(gè)平表面,在傳統(tǒng)法中利用L方法[11]進(jìn)行多視圖的重建,L方法可以獲得全局最優(yōu)解,因此將這種方法記為“L”法。實(shí)驗(yàn)中采用的性能參數(shù)為:重投影均方根誤差為i=1,2,…,M,j=1,2,…,N,xij代表檢測(cè)點(diǎn)坐標(biāo),代表反投影點(diǎn)坐標(biāo),單位為像素;空間位置均方根誤差為
仿真實(shí)驗(yàn)中需要知道產(chǎn)生檢測(cè)數(shù)據(jù)的真正平表面參數(shù),設(shè)置如下:將平表面定為Xw-Yw平表面,即Zw=0平表面,隨機(jī)生成范圍在[-1,1]內(nèi)的100個(gè)點(diǎn),即N=100,要保證這些點(diǎn)在每個(gè)攝像機(jī)下都是可視的。假設(shè)用6個(gè)攝像機(jī)進(jìn)行拍攝,實(shí)驗(yàn)中每臺(tái)攝像機(jī)的內(nèi)參數(shù)設(shè)為
設(shè)定每個(gè)攝像機(jī)之間大概相差10°,圖像分辨率為640×480。仿真實(shí)驗(yàn)環(huán)境如圖4所示,圖4中用圓圈來表示攝像機(jī)的位置,點(diǎn)云為隨機(jī)生成的Zw=0平表面上的點(diǎn)。
圖4 仿真實(shí)驗(yàn)環(huán)境示意圖Fig.4 The schematic diagram of simulation environment
在每個(gè)圖像點(diǎn)上增加零均值、標(biāo)準(zhǔn)差為σ的高斯噪聲,噪聲大小分別為0.2、0.4、0.6、0.8、1.0個(gè)像素,對(duì)每種噪聲分別進(jìn)行100次實(shí)驗(yàn)。由于對(duì)帶有噪聲的數(shù)據(jù)都進(jìn)行了100次實(shí)驗(yàn),每次實(shí)驗(yàn)都會(huì)產(chǎn)生一個(gè)重建結(jié)果,數(shù)據(jù)過多,不再列出,只給出3種平表面重建方法的性能比較,如表1所示。從表1可以看出,RPE-GA法和TE-GA法的精度基本一致,可以明顯看出空間位置均方根誤差比較小,這是因?yàn)槭峭ㄟ^最小化的空間點(diǎn)位置誤差進(jìn)行優(yōu)化求解的,所以空間位置均方根誤差較?。欢鳯方法的重投影均方根誤差較小,這是因?yàn)樵贚法中是通過最小化圖像平面上的重投影誤差進(jìn)行優(yōu)化求解的,所以重投影均方根誤差較小。由于在L法中沒有考慮到景物平表面的約束,所以空間位置均方根誤差較大,也就與真實(shí)的景物平表面相差較遠(yuǎn)。
表1 3種平表面重建方法的性能比較
真實(shí)圖像實(shí)驗(yàn)采用Zhang從不同角度拍攝平面模版的5幅圖像(http://research.microsoft.com/en-us/um/people/zhang/calib/),如圖5所示,圖像分辨率大小為640×480,平面模版上共有256個(gè)角點(diǎn),即N=256。
圖5 平面模版的5幅圖像Fig.5 Five images of a planer pattern
真實(shí)圖像的重建結(jié)果如圖6所示。從圖6可以看出本文方法的重建結(jié)果基本上位于一個(gè)平表面上,而傳統(tǒng)的L法重建結(jié)果非常散亂,并不位于一個(gè)平表面上。真實(shí)圖像實(shí)驗(yàn)中3種重建平表面方法的性能比較如表2所示。同樣的,RPE-GA法和TE-GA法的空間位置均方根誤差比較小,與真實(shí)的景物平表面相接近;而L法的重投影均方根誤差較小,但空間位置均方根誤差卻比較大,也就是與真實(shí)的景物平表面相差較遠(yuǎn)。從原理上來講,RPE-GA和TE-GA 2種優(yōu)化方法是一樣的,都是將空間點(diǎn)限制在一個(gè)平表面上,然后最小化空間點(diǎn)之間的距離誤差,所以它們的精度基本一致。只是由于RPE-GA法中多了參數(shù)αij,所以RPE-GA法的計(jì)算會(huì)比TE-GA法復(fù)雜。本文方法的優(yōu)越性:1)利用三維空間點(diǎn)都位于同一個(gè)景物平表面這一先驗(yàn)知識(shí),計(jì)算出的平表面精度較高;2)可以推廣到多幅視圖。
需要指出的是RPE-GA法和TE-GA法對(duì)局外點(diǎn)都比較敏感,事實(shí)上,目前已有的基于兩視圖景物平表面重建方法,如Olsson[7]提出的重建方法對(duì)局外點(diǎn)也比較敏感。如果用在場景比較復(fù)雜的真實(shí)圖像實(shí)驗(yàn)中,會(huì)存在2類局外點(diǎn),一類是在匹配錯(cuò)誤的點(diǎn),一類是不位于一個(gè)平表面上的點(diǎn)。首先需要采用一些剔除局外點(diǎn)的方法對(duì)局外點(diǎn)進(jìn)行剔除,Olsson采用L1方法對(duì)局外點(diǎn)進(jìn)行了剔除,但是從實(shí)驗(yàn)結(jié)果中可以很明顯的看出,即使對(duì)局外點(diǎn)進(jìn)行了剔除,仍然有很多圖像點(diǎn)不位于同一個(gè)平表面上,用這些不位于同一個(gè)平表面上的圖像點(diǎn)去重建一個(gè)平表面勢(shì)必是不準(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
在傳統(tǒng)的三維景物中平表面重建法中,都是先重建出三維點(diǎn),再用這些三維點(diǎn)去擬合景物平表面,由于沒有利用三維空間點(diǎn)位于一個(gè)平表面的先驗(yàn)知識(shí),所以重建出的三維點(diǎn)并不位于一個(gè)平表面上。針對(duì)這個(gè)問題,本文提出了2種三維景物中平表面重建模型,2種模型都是考慮到重建的空間點(diǎn)應(yīng)位于同一個(gè)平面內(nèi),以此為約束條件分別在空間平面上進(jìn)行最小化反投影誤差和最小化轉(zhuǎn)移誤差,再采用GA算法對(duì)2種模型進(jìn)行了優(yōu)化求解,從而獲得平表面重建結(jié)果。實(shí)際上,基于最小化反投影誤差的平表面重建法和基于最小化轉(zhuǎn)移誤差法的平表面重建法基本原理相同,只是計(jì)算復(fù)雜度不同。實(shí)驗(yàn)表明了本文方法的優(yōu)越性。另外如何剔除掉局外點(diǎn),或者說如何從復(fù)雜場景中提取出位于一個(gè)平表面上的圖像點(diǎn),將RPE-GA法和TE-GA法用于復(fù)雜場景的真實(shí)圖像實(shí)驗(yàn)上,將是未來需要研究的問題。據(jù)知,目前在國內(nèi)外考慮三維點(diǎn)位于空間景物平面的約束條件,針對(duì)多視圖的三維景物中平表面重建的研究還處于空白。
參考文獻(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.