門朝光，邊繼龍，李香

(哈爾濱工程大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，黑龍江 哈爾濱，150001)

立體匹配是計(jì)算機(jī)視覺領(lǐng)域中的一個(gè)熱點(diǎn)問(wèn)題，它通過(guò)匹配2幅或多幅圖像獲得對(duì)應(yīng)點(diǎn)視差，再通過(guò)三角測(cè)量原理計(jì)算出景物的深度，在計(jì)算機(jī)視覺、立體測(cè)繪、航空攝影測(cè)量、遙感等領(lǐng)域有著廣泛的應(yīng)用。目前，雖然基于局部和基于全局的立體匹配方法[1]取得了一定進(jìn)展，但由于匹配過(guò)程是一個(gè)病態(tài)過(guò)程，并且存在著大量的不確定因素，使得立體匹配仍然是一個(gè)活躍的研究領(lǐng)域。由于重復(fù)紋理、非紋理區(qū)域和遮擋的存在使像素點(diǎn)之間的匹配具有很大程度的不確定性，特別當(dāng)立體像對(duì)不是在同一時(shí)間獲取時(shí)，移動(dòng)或消失的目標(biāo)會(huì)使匹配變得更加困難。為解決上述因素的影響，在立體視覺領(lǐng)域中出現(xiàn)了另一個(gè)分支即小基高比立體匹配[2-8]。基于小基高比的立體匹配是指在攝影基線與攝像機(jī)高度比值較小，對(duì)目標(biāo)物體進(jìn)行拍照獲得立體像對(duì)，再利用立體匹配技術(shù)獲得對(duì)應(yīng)點(diǎn)視差，最后，根據(jù)三角測(cè)量原理計(jì)算出物體的高程信息?；谛』弑鹊牧Ⅲw匹配算法中存在 1個(gè)相關(guān)基本等式[2-4](Central equation of correlation)，此等式把測(cè)量視差和真實(shí)視差聯(lián)系起來(lái)，可在匹配結(jié)束后對(duì)相關(guān)基本等式進(jìn)行求解以恢復(fù)其真實(shí)視差。通過(guò)求解此等式可以解決小基高比立體匹配中存在的邊界膨脹問(wèn)題即黏合現(xiàn)象。相關(guān)基本等式求解問(wèn)題屬于數(shù)學(xué)上的反問(wèn)題，具有不適定性[9]，需要對(duì)相關(guān)基本等式進(jìn)行正則化處理，使問(wèn)題變得適定可解。本文通過(guò)在目標(biāo)泛函中引入全變差項(xiàng)進(jìn)行正則化處理，把等式求解問(wèn)題轉(zhuǎn)化為求解目標(biāo)泛函的極值問(wèn)題。目前，數(shù)字圖像處理領(lǐng)域中變分法受到了許多研究者的關(guān)注，它廣泛應(yīng)用于圖像復(fù)原[10-11]、圖像融合[12]、圖像分割[13]和立體匹配[14]中?；谧兎址ǖ膬?yōu)勢(shì)在于它能提供一個(gè)良好的數(shù)學(xué)形式體系，而且可以生成亞像素級(jí)視差。變分法中的正則化參數(shù)選取關(guān)系到解的逼近性與數(shù)值穩(wěn)定性，其選取方法可分為先驗(yàn)和后驗(yàn)2類策略。其中，先驗(yàn)策略是在求出正則解之前就已將正則參數(shù)確定。這種方法往往具有理論分析價(jià)值，在實(shí)際中常常難以驗(yàn)證其賴以施用的條件；后驗(yàn)策略則是在計(jì)算正則解的過(guò)程中根據(jù)一定的原則確定與原始數(shù)據(jù)的誤差水平相匹配的正則參數(shù)，此方法比較實(shí)用。本文采用后驗(yàn)策略根據(jù)Morozov偏差原理推導(dǎo)出正則參數(shù)選擇公式，并設(shè)計(jì)一個(gè)迭代過(guò)程進(jìn)行求解。根據(jù)Morozov原理求解正則參數(shù)需要了解系統(tǒng)的誤差水平，然而，在立體匹配中我們無(wú)法事先獲知該誤差水平。若根據(jù)經(jīng)驗(yàn)設(shè)置此誤差水平，則在匹配不同的立體像對(duì)時(shí)需要選擇不同的經(jīng)驗(yàn)值，這會(huì)影響算法的適應(yīng)性，為此本文根據(jù)啟發(fā)式信息估計(jì)該誤差水平。最后，根據(jù)延遲擴(kuò)散定點(diǎn)迭代思想[15](lagged diffusivity fixed point iteration，F(xiàn)P)獲得目標(biāo)泛函的迭代傳播等式，通過(guò)共軛梯度法進(jìn)行求解。

1 相關(guān)性分析

假設(shè)u和u?代表立體像對(duì)中左右圖像，并滿足經(jīng)典立體模型：

其中：函數(shù) λ(x)在空間上緩慢變化；視差函數(shù) ε(x)是一有界函數(shù)，描述立體像對(duì)u和之間的幾何畸變。在獲取立體像對(duì)的過(guò)程中，若2次成像之間的角度較大即大基高比，則立體像對(duì)會(huì)受光照變化(例如云層遮擋)、物體遮擋和非朗伯平面等因素的影響而導(dǎo)致它們之間存在較大的輻射差異，致使假設(shè)立體模型存在較大的偏差；若2次成像之間的角度較小即小基高比，則可以認(rèn)為2次成像時(shí)的光照條件近似一樣而且小角度可以減少物體遮擋，減弱非朗伯平面對(duì)輻射差異的

此立體模型假設(shè)立體像對(duì)u和u?之間僅存在幾何畸變，不存在輻射差異，這種模型只有在基高比 b/h較小時(shí)才能成立，而且基高比越小模型越精確。

基于小基高比的立體匹配方法[2-8]屬于局部匹配方法，該方法的匹配精度受支撐窗口大小影響。當(dāng)窗口大小增加時(shí)匹配精度隨之增加，但同時(shí)大窗口也會(huì)導(dǎo)致算法在視差不連續(xù)處產(chǎn)生誤匹配，在立體匹配中把這種誤匹配稱為黏合現(xiàn)象，這種現(xiàn)象會(huì)導(dǎo)致視差圖中的前景區(qū)膨脹失真。如圖1所示，左圖像中存在3個(gè)點(diǎn)a，b和c，它們分別與右圖像中的′，′和′ 3個(gè)點(diǎn)相對(duì)應(yīng)，其視差分別為0，db和dc。在匹配時(shí)由于b點(diǎn)的支撐窗口包含了建筑物邊界致使算法在b點(diǎn)產(chǎn)生錯(cuò)誤的視差db=dc，從而導(dǎo)致前景區(qū)的建筑物產(chǎn)生膨脹。影響，在這種條件下生成的立體像對(duì)可以近似滿足這一模型。在小基高比條件下，立體模型可以簡(jiǎn)化為：

圖1 黏合現(xiàn)象Fig.1 Adhesion phenomenon

假設(shè)在小基高比模型(2)中真實(shí)視差 ε(x)足夠小，通過(guò)測(cè)量視差對(duì)這個(gè)模型進(jìn)行一階近似，可以發(fā)現(xiàn)通過(guò)相關(guān)基本等式(4)能將最大化式(3)所獲得的測(cè)量視差與真實(shí)視差聯(lián)系起來(lái)。

式中：函數(shù) φ ( x0- x )表示以x0為中心的窗口函數(shù)；ε(x)表示真實(shí)視差函數(shù)；0()d x為圖像u在x0點(diǎn)的相關(guān)密度函數(shù)[4]，其數(shù)學(xué)表達(dá)式為：

相關(guān)密度函數(shù)只依賴于參考圖像 u，它表達(dá)了窗口中心x0周圍的邊緣方向。通過(guò)求解相關(guān)基本等式可以解決小基高比立體匹配中存在的黏合現(xiàn)象。Delon等[4]通過(guò) delta函數(shù)近似式(4)中的相關(guān)密度函數(shù)將等式(4)簡(jiǎn)化為ε(x1)= m(x0)，再通過(guò)質(zhì)心校正公式(6)計(jì)算出質(zhì)心的灰度：

然后，在局部窗口內(nèi)尋找一個(gè)與此灰度x1相似的像素位置作為質(zhì)心位置并將測(cè)量視差賦予此位置。雖然通過(guò)簡(jiǎn)單的質(zhì)心校正取得了一定的成效，但根據(jù)該校正公式計(jì)算的灰度確定質(zhì)心位置會(huì)導(dǎo)致誤差傳播。為此，本文提出通過(guò)變分原理求解該等式來(lái)解決小基高比立體匹配中存在的黏合現(xiàn)象，該方法可以有效減少匹配中存在的黏合現(xiàn)象，提高視差圖的準(zhǔn)確率。

2 模型分析及求解

2.1 模型分析

小基高比立體匹配中存在的相關(guān)基本等式(4)可以表示成形如式(7)的第一類算子方程，其數(shù)學(xué)表達(dá)式為：

由于式(7)中的算子K是1個(gè)緊算子，因此，導(dǎo)致求解相關(guān)基本等式問(wèn)題成為1個(gè)病態(tài)問(wèn)題。為使等式變得適定可解，經(jīng)常使用Tikhonov正則化方法對(duì)等式進(jìn)行正則化處理。根據(jù)Tikhonov正則化方法，將求解等式(7)問(wèn)題轉(zhuǎn)化為求解目標(biāo)泛函極值問(wèn)題，其目標(biāo)泛函為：

式中：?ε為ε的梯度，其數(shù)學(xué)表達(dá)式為

圖2 懲罰函數(shù)的比較Fig.2 Comparison of different penalty functions

第1種懲罰函數(shù)式(9)(如圖2(a)所示)在視差函數(shù)ε不連續(xù)時(shí)經(jīng)常導(dǎo)致虛假振蕩(spurious oscillations)；第2種懲罰函數(shù)式(10)要求視差函數(shù)ε處處平滑，然而在物體邊界(即視差不連續(xù)處)常常會(huì)違背這一假設(shè)，因此，這2種懲罰函數(shù)都不太適合此應(yīng)用。為保持函數(shù)的不連續(xù)性，防止在邊緣處過(guò)度平滑，許多研究者引入了全變差[10-12]來(lái)解決此問(wèn)題，其全變差的數(shù)學(xué)表達(dá)式為：

此懲罰函數(shù)的形狀如圖2(b)所示。由于此懲罰函數(shù)不要求視差函數(shù)ε處處平滑，因此該函數(shù)是一個(gè)很好的候選正則函數(shù)。然而該函數(shù)的歐式范數(shù)在0處不可微，為避免這一問(wèn)題，正則函數(shù)修改為：

其修改后的懲罰函數(shù)形狀如圖 2(c)所示。在 β=0時(shí)，可以簡(jiǎn)化為全變差函數(shù)

根據(jù)修改的全變差懲罰函數(shù)式(13)相應(yīng)的目標(biāo)泛函可以表達(dá)為：

然后，通過(guò)對(duì)目標(biāo)泛函E(ε)求導(dǎo)獲得該泛函的歐拉-拉格朗日方程：

式中：K*代表伴隨算子，其數(shù)學(xué)表達(dá)式為

L(ε)為微分算子，它對(duì)函數(shù)ω的作用可以通過(guò)下式表示：

最后，根據(jù)Vogel提出的延遲擴(kuò)散定點(diǎn)迭代思想最小化目標(biāo)泛函式(14)。其FP迭代可以表達(dá)為：

每次迭代必須通過(guò)求解線性傳播式(18)來(lái)獲得下一步迭代εk+1，迭代過(guò)程的傳播率取決于前一迭代εk。

2.2 數(shù)值求解

經(jīng)過(guò)離散化處理之后，式(18)轉(zhuǎn)化為超大稀疏線性方程組：

由于系數(shù)矩陣A為超大稀疏矩陣，直接通過(guò)高斯消除法求解線性方程組的方法不可行，因此，本文通過(guò)迭代的共軛梯度法來(lái)求解線性方程組，其共軛梯度法如下。

算法1(共軛梯度法[16])：

k=0

r0=b-Ax0

while rk≠0

k=k+1

if k=1

p1=r0

else

βk=(rk-1)Trk-1/(rk-2)Trk-2

pk= rk-1+βkpk

end

αk=(rk-1)Trk-1/(pk)TApk

xk=xk-1+αkpk

rk=rk-1-αkApk

end

x=xk

3 正則化處理

3.1 正則參數(shù)選擇

由于正則化參數(shù)α的選取關(guān)系到解的逼近性與數(shù)值穩(wěn)定性，為此，根據(jù)Morozov原理[17]選擇此參數(shù)以獲得較為精確的視差。Morozov原理闡述了怎樣利用觀測(cè)數(shù)據(jù)選擇正則參數(shù) α，它在線性反問(wèn)題領(lǐng)域中受到了很大關(guān)注。由于正則化產(chǎn)生的誤差等于觀測(cè)數(shù)據(jù)所引起的誤差，因此，可以根據(jù)這一原理選擇正則參數(shù)α：

在計(jì)算正則參數(shù)的過(guò)程中需要利用以下幾個(gè)法則，本文在此給出這些法則和其證明過(guò)程。令F(α)是關(guān)于正則參數(shù)的函數(shù)，其數(shù)學(xué)表達(dá)式為：

法則1 對(duì)任意α＞0，式(21)存在惟一最小化解，它可以描述成下面系統(tǒng)的解：

或者用內(nèi)積可表示為：

法則2 F(α)的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)可以表示為：

證明：

根據(jù)式(23)可得：

因此，可得：

法則 3 F(α)是關(guān)于正則參數(shù)的函數(shù)，對(duì)任意α＞0，存在如下等式：

證明 根據(jù)法則1，有如下等式成立：

對(duì)上述等式兩邊關(guān)于α求導(dǎo)得：

根據(jù)法則2，上述等式可以轉(zhuǎn)化為：

對(duì)等式兩邊積分得：

根據(jù) F(α)和 F′(α)可以將式(20)轉(zhuǎn)化為：

為了獲得合理的正則參數(shù)，根據(jù)Morozov原理推導(dǎo)出正則參數(shù)選擇方法。為推導(dǎo)模型函數(shù)，對(duì)式(25)的第3項(xiàng)進(jìn)行如下近似：

式中：T1是待定常數(shù)。通過(guò)上述近似，式(25)可以簡(jiǎn)化為：

對(duì)等式(28)兩邊積分可得：

為簡(jiǎn)化正則參數(shù)α的計(jì)算，令m(0)=0，其模型函數(shù)m(α)的表達(dá)式可以簡(jiǎn)化為：

對(duì)模型函數(shù)m(α)關(guān)于α求導(dǎo)得：

最后，通過(guò)迭代方法求解正則參數(shù)α，具體算法如下。

(1) 根據(jù)式(21)和式(24)計(jì)算 F(αk)和 F′(αk)，再根據(jù)式(32)計(jì)算Ck和Tk。

通過(guò)求解式(32)可得：

(2) 根據(jù)式(30)和式(31)計(jì)算 m(αk)與 m'(αk)。

(3) 根據(jù) Morozov偏差方程(26)計(jì)算正則參數(shù)αk+1。

3.2 估計(jì)誤差水平

計(jì)算正則參數(shù)時(shí)需要使用觀測(cè)誤差水平δ，然而，立體匹配中往往無(wú)法事先獲知該誤差水平。為了能精確的估計(jì)正則參數(shù)，在計(jì)算正則參數(shù)之前首先利用啟發(fā)式方法估計(jì)該誤差水平。根據(jù)Morozov原理有如下不等式關(guān)系成立：

利用式(34)可推導(dǎo)出誤差觀測(cè)水平的存在區(qū)間：

在數(shù)值實(shí)現(xiàn)過(guò)程中，利用模型函數(shù)m(α)預(yù)測(cè)觀測(cè)誤差δ：

具體實(shí)現(xiàn)過(guò)程如下：

(1) 給定比例常數(shù) σ ∈ [ 0 .5,1]和初始正則參數(shù) α0＞0。

(2) 利用式(21)和式(24)計(jì)算 F(α0)和 F'(α0)。然后，利用式(37)計(jì)算式(30)中的參數(shù)C和T。

再利用式(38)計(jì)算正則參數(shù)α1：

式中：y0為曲線 y=m(α)上點(diǎn)(α0,m(α0))的切線與 y 軸的截距，其數(shù)學(xué)表達(dá)式為

(3) 根據(jù)步驟(2)中計(jì)算的正則參數(shù)α1，利用式(21)和式(24)計(jì)算 F(α1)和 F'(α1)，再利用式(25)計(jì)算參數(shù)C0，然后，利用式(29)和式(40)計(jì)算參數(shù)C1和T1。

利用式(29)計(jì)算m(0)和m(1)，誤差δ表示為：

3.3 算法總體流程

基于迭代傳播的小基高比立體匹配方法是在Delon等[4]提出的基于多分辨率的小基高比立體匹配方法的基礎(chǔ)上通過(guò)對(duì)匹配中存在的相關(guān)基本等式進(jìn)行迭代求解來(lái)恢復(fù)其真實(shí)視差值。該方法可以減少小基高比立體匹配過(guò)程中因匹配窗口跨越邊界而導(dǎo)致的黏合現(xiàn)象，立體匹配中的這種黏合現(xiàn)象會(huì)導(dǎo)致視差圖中前景區(qū)產(chǎn)生“膨脹”。算法具體流程如下：

(1) 利用小基高比立體匹配算法[4]獲得初始測(cè)量視差ε0，并給定比例常數(shù)σ和初始正則參數(shù)α0。

(2) 利用3.2節(jié)中的算法獲得誤差水平δk。

(3) 利用3.1節(jié)中的算法獲得正則參數(shù)αk。

(4) 利用共軛梯度法求解線性方程組(19)獲得視差 εk+1，若或者 k ≥kmax則停止迭代；否則，k=k+1，轉(zhuǎn)入步驟(2)繼續(xù)迭代。

4 實(shí)驗(yàn)分析

本實(shí)驗(yàn)使用 CNES提供的航空攝影像對(duì)Toulouse(如圖3(a)，(b)所示)對(duì)本文算法的合理性進(jìn)行驗(yàn)證。這幅立體像對(duì)是1幅航空影像圖，基高比約為0.045，地面分辨率為 R=0.5，其獲取時(shí)間間隔為 20 min，這導(dǎo)致立體像對(duì)中存在明顯的運(yùn)動(dòng)和陰影移動(dòng)，增加了視差估計(jì)的難度。實(shí)際上，理想的小基高比立體像對(duì)的2次獲取時(shí)間間隔較小，這樣可以有效的減少立體像對(duì)間的輻射差異和景物運(yùn)動(dòng)對(duì)匹配產(chǎn)生的不利影響。

圖3 Toulouse立體像對(duì)Fig.3 Toulouse stereo images

理論上講，小基高比條件下形成的立體像對(duì)非常相似，具有較小的視差搜索范圍。這幅立體像對(duì)的視差范圍為[-2，2]，可近似模擬小基高比條件下的立體像對(duì)，這種相似性使匹配過(guò)程變得更容易。圖 3(a)和(b)所示為這幅立體像對(duì)的左右圖像，看出它們非常相似，具有較小的幾何畸變與視差搜索范圍。其真實(shí)視差圖如3(c)所示。圖4(a)為Delon等[4]提出的基于多分辨率的小基高比立體匹配方法獲得的實(shí)驗(yàn)結(jié)果，通過(guò)與真實(shí)視差圖對(duì)比可以看出該方法已經(jīng)獲得了較好的效果，但是，在背景處還是存在一些明顯的誤匹配，而且在物體邊界處伴隨著較多黏合現(xiàn)象(如圖 4(b)所示)。圖4(c)所示為由本文算法生成的視差圖。從圖4可見：由本文算法生成的視差圖中場(chǎng)景各主要組成部分都清晰地檢測(cè)出來(lái)，特別是建筑物四周只發(fā)生微小變化的那些部分也都被清晰地檢測(cè)出來(lái)。通過(guò)對(duì)比圖4(b)和 4(d)可以看出：本文算法的黏合現(xiàn)象要明顯少于文獻(xiàn)[4]中的黏合現(xiàn)象，而且該算法生成的視差圖的準(zhǔn)確率高達(dá)95%以上。

為驗(yàn)證小基高比立體匹配的高程精度，本實(shí)驗(yàn)使用由空間機(jī)電研究所提供的1幅帶有真實(shí)高程信息的小基高比立體像對(duì)(如圖 5(a)和(b)所示)。該立體像對(duì)的基高比約為0.05，地面分辨率為0.2。在該實(shí)驗(yàn)中，根據(jù)本文算法生成的視差圖(如圖 5(c)所示)計(jì)算幾個(gè)主要建筑物的高程信息，這些建筑被標(biāo)記在圖5(d)中，對(duì)應(yīng)建筑物的視差、計(jì)算高程和真實(shí)高程如表1所示。高程計(jì)算公式如下：式中：δ為視差；B/H為基高比；R為地面分辨率。從表1可知：該算法具有較高的高程精度，其平均高程誤差為0.33 m，像元匹配差異精度為0.082 5，優(yōu)于1/10個(gè)像元。

圖4 黏合現(xiàn)象比較Fig.4 Comparison of adhesion phenomenon

表1 高程信息Table 1 Information of DEM

圖5 高程精度驗(yàn)證Fig.5 Verification of DEM accuracy

5 結(jié)論

(1) 根據(jù)啟發(fā)式信息估計(jì)視差誤差水平，解決了先驗(yàn)參數(shù)估計(jì)問(wèn)題，使匹配算法具有更強(qiáng)的適應(yīng)性。

(2) 所提出的迭代正則參數(shù)選擇方法解決了正則解的穩(wěn)定性與逼近性問(wèn)題，使匹配算法獲得了更加精確的視差。

(3) 通過(guò)共軛梯度法求解目標(biāo)泛函的迭代傳播等式加快了匹配算法的收斂速度。

(4) 該方法能有效地減少小基高比立體匹配中的黏合現(xiàn)象，提高視差圖的準(zhǔn)確率，其視差圖的準(zhǔn)確率可達(dá)95%以上，且亞像素精度可優(yōu)于1/10個(gè)像元。

[1] 羅三定, 陳海波. 基于區(qū)域增長(zhǎng)的自適應(yīng)窗口立體匹配算法[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 2005, 36(6): 1042-1047.LUO San-ding, CHEN Hai-bo. Stereo matching algorithm of adaptive window based on region growing[J]. Journal of Central South University: Science and Technology, 2005, 36(6):1042-1047.

[2] Delon J. Fine comparison of images and other problems[D].Paris: ENS Cachan. Mathematics Department, 2004: 89-108.

[3] Facciolo G. Variational adhesion correction with image based regularization for digital elevation models[D]. Montevicleo:University of the Republic. Engineering College, 2005: 37-49.

[4] Delon J, Rouge B. Small baseline stereovision[J]. Journal of Mathematical Imaging and Vision, 2007, 28(3): 209-223.

[5] Igual L, Preciozzi J, Garrido L, et al. Automatic low baseline stereo in urban areas[J]. Inverse Problems and Imaging, 2007,1(2): 319-348.

[6] Morgan G L K, LIU Jian-guo, YAN Hong-shi. Sub-pixel stereo matching for DEM generation from narrow baseline stereo imagery[C]//Proceedings of International Geoscience and Remote Sensing Symposium. Boston: IEEE, 2008: 1284-1287.

[7] Sabater N, Blanchet G, Moisan L,et al. Review of low-baseline stereo algorithms and benchmarks[C]//Proceedings of Image and Signal Processing for Remote Sensing XVI. Toulouse: IEEE,2010: 1-12.

[8] Morgan G L K, LIU Jian-guo, YAN Hong-shi. Precise subpixel disparity measurement from very narrow baseline stereo[J].IEEE Transactions on Geoscience and Remote Sensing, 2010,48(9): 3424-3433.

[9] 王彥飛. 反問(wèn)題的計(jì)算方法及其應(yīng)用[M]. 北京: 高等教育出版社, 2007: 3-15.WANG Yan-fei. Computational method for inverse problem and their application[M]. Beijing: High Education Press, 2007: 3-15.

[10] Aujol J F. Some first-order algorithms for total variation based image restoration[J]. Mathematical Imaging and Vision, 2009,34(3): 307-327.

[11] Beck A, Teboulle M. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[J]. IEEE Transactions on Image Processing, 2009,18(11): 2419-2434.

[12] Mrityunjay K, Sarat D. A total variation-based algorithm for pixel-Level image fusion[J]. IEEE Transactions on Image Processing, 2009, 18(9): 2137-2143.

[13] Chung G, Vese L A. Image segmentation using a multilayer level-set approach[J]. Computing and Visualization in Science,2009, 12(6): 267-285.

[14] Benari R, Sochen N. Stereo matching with Mumford-Shah regularization and occlusion handling[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2010, 32(11):2071-2084.

[15] Vogel C R, Oman M E. Iterative methods for total variation denoising[J]. SIAM Journal on Scientific Computing, 1996,17(1): 227-238.

[16] Golub G H, Loan C F V. Matrix computations[M]. London:Johns Hopkins Univ, 1996: 520-532.

[17] Kunisch K, ZOU Jun. Iterative choices of regularization parameters in linear inverse problems[J]. Inverse Problems, 1998,14: 1247-1264.