伍秋玉，張明新，劉永俊，鄭金龍

(1.中國礦業(yè)大學(xué) 計算機科學(xué)與技術(shù)學(xué)院，江蘇 徐州 221116； 2.常熟理工學(xué)院 計算機科學(xué)與工程學(xué)院，江蘇 常熟 215500；3.東北大學(xué) 計算機科學(xué)與工程學(xué)院，沈陽 110819)(*通信作者電子郵箱zmxcslg@163.com)

0 引言

在微納米計算機視覺中，基于視覺的微納米尺度3D形貌重建，對于較全面地理解樣品特性和對操作過程進行評估具有重要意義。比較常用的微觀3D重建方法[1]主要有體積恢復(fù)方法、立體視覺深度恢復(fù)方法、聚焦深度恢復(fù)方法以及離焦深度恢復(fù)方法。與其他的3D重建方法相比，離焦深度恢復(fù)方法由于所需圖片少、設(shè)備簡單、操作便利等優(yōu)點，近些年得到廣泛關(guān)注與深入研究[2]。離焦深度恢復(fù)最早由Pentland[3]提出，該方法利用二維圖像的離焦程度特征與景物深度之間的映射關(guān)系，反解出景物的三維深度信息[4-5]。

在景物的離焦深度恢復(fù)方法中，首先必須獲取景物不同程度的離焦圖像，這就導(dǎo)致要改變攝像機參數(shù)。但是，在微納米圖像觀測中，觀測空間非常有限，而且使用的是具有高放大倍數(shù)的相機，所以相機的成像模型會隨著攝像機參數(shù)的改變而改變[6]。受限于微納米觀測中的條件，魏陽杰等[7]提出了一種新的基于參數(shù)固定的單目視覺攝像機的三維離焦形貌恢復(fù)方法。該方法利用相對模糊度[7-8]及熱輻射方程[7,9]來求解景物的深度信息，并將深度信息的求解轉(zhuǎn)化成動態(tài)優(yōu)化問題來求解。

在求解離焦深度恢復(fù)動態(tài)優(yōu)化問題時，文獻[7]采用經(jīng)典的迭代收縮閾值算法(Iterative Shrinkage Thresholding Algorithm， ISTA)[10]來求解，但ISTA是梯度下降法的延伸，迭代過程僅考慮當(dāng)前點的信息進行梯度估計更新迭代點，優(yōu)化過程呈“之”字形[10-12]向極小值點靠近，收斂速度比較慢；而且該算法在迭代過程中采用固定步長，在靠近極小值點時，此時的固定步長可能大于實際的迭代步長，導(dǎo)致算法迭代效率不佳，從而使得重建的微觀3D形貌精度不高。

針對ISTA求解離焦深度恢復(fù)動態(tài)優(yōu)化問題的不足，本文提出了一種基于加速算子梯度估計和割線線性搜索的方法優(yōu)化ISTA——FL-ISTA(Fast Linear ISTA)。該算法在ISTA的基礎(chǔ)上加入一個加速算子，該加速算子由當(dāng)前點和前一個點的線性組合構(gòu)成，利用加速算子重新進行梯度估計，更新迭代點，以加快算法的迭代速度；而且為了改變ISTA迭代步長固定的限制，結(jié)合割線線性搜索尋找最優(yōu)迭代步長，以提高算法的迭代效率。將FL-ISTA應(yīng)用于標準500 nm尺度柵格的深度信息重建，實驗結(jié)果表明，與ISTA、快速迭代收縮閾值算法(Fast ISTA, FISTA)和單調(diào)快速迭代收縮閾值算法(Monotohy FISTA, MFISTA)相比，F(xiàn)L-ISTA收斂速度更快，重建的深度信息值下降更快，更接近標準500 nm柵格尺度，并且相對誤差、平均誤差、均方差更小，有效地提高了求解離焦深度恢復(fù)動態(tài)優(yōu)化問題的效率和重建的微觀3D形貌的精度。

1 全局離焦深度恢復(fù)

1.1 離焦成像模型

在圖像處理中，圖像的模糊化通常可以用卷積的形式表示如下：

I2(x,y)=I1(x,y)*H(x,y)

(1)

其中:I1(x,y)、I2(x,y)和H(x,y)分別代表清晰圖像與模糊圖像以及點擴散函數(shù)，*表示卷積。根據(jù)點擴散原理，點擴散函數(shù)可由二維高斯函數(shù)來近似表示，并且其中的模糊擴散參數(shù)σ表示圖像的模糊程度值。由于熱輻射方程具有各向同性，所以式(1)又可以表示為：

(2)

并且:

σ2=2tc

(3)

如果深度圖是理想平面，那么c是一個常數(shù)，否則可以表示為如下形式:

(4)

其中▽和▽·分別代表梯度算子和微分算子：

由上式可知，首先需要已知清晰圖像g(x,y)來求解熱輻射方程，但這是一個復(fù)雜的過程，因此，F(xiàn)avaro等[8]提出了相對模糊的概念。

假設(shè)有兩張不同程度的模糊圖像I1(x,y)和I2(x,y)，模糊擴散系數(shù)分別為σ1和σ2，并且σ1<σ2，那么I2(x,y)可以寫成如下形式:

(5)

(6)

而式(4)又可以寫成如下形式:

(7)

并且在Δt時刻，u(x,y,t)=I2(x,y)，Δt可被定義為:

Δσ2=2(t2-t1)c?2Δtc

(8)

而且：

(9)

其中:ηi(i=1,2)表示模糊圓半徑，而λ則表示模糊度與模糊圓半徑之間的常數(shù)。

(10)

其中:v、f、s分別為攝像機的像距、焦距和物距，D為凸透鏡半徑。

1.2 深度信息重建

假設(shè)圖像I1(x,y)是物距變化之前的離焦圖像，它的深度信息為S1(x,y)，而圖像I2(x,y)是物距變化之后的離焦圖像，它的深度信息為S2(x,y)，已知理想焦距為s0，在這一部分，圖像I2(x,y)的深度信息是由物距變化ΔS來獲得，原理示意圖如圖1所示。

首先由以上條件建立熱輻射方程:

(11)

其中相對模糊度又可以表示為:

(12)

圖1 深度信息原理圖

本文定義:

(13)

因此，最終的深度信息為:

(14)

為了更好地解決離焦深度恢復(fù)動態(tài)優(yōu)化問題，建立如下目標函數(shù)來計算熱輻射方程:

(15)

由于上述的求解過程可能是病態(tài)的，所以在目標函數(shù)最后加上一個Tikhonov正則項，表示為:

ρ‖▽s2(x,y)‖2+ρl‖s2(x,y)‖2

(16)

上式第三項是關(guān)于深度圖的平滑約束，第四項是保證深度圖的有界性，實際上能量系數(shù)ρ>0，l>0。損失的能量表示為:

E(s)=?(z(x,y,Δt)-I2(x,y))2dxdy+

ρ‖▽s‖2+ρl‖s‖2

(17)

因此，可以通過求解如下優(yōu)化問題來求解景物深度信息:

(18)

s.t. (11)和(14)

2 FL-ISTA

2.1 ISTA求解離焦深度恢復(fù)動態(tài)優(yōu)化問題

對于1.2節(jié)中的離焦深度恢復(fù)動態(tài)優(yōu)化問題，采用ISTA進行求解，令

E1(s)=?(z(x,y,Δt)-I2(x,y))2dxdy

(19)

于是，對于無約束最優(yōu)化問題:

E1(s)=?(z(x,y,Δt)-I2(x,y))2dxdy

解決上式的梯度算法為:

sj=sj-1-tj▽E1(sj-1)

(20)

其中tj表示迭代步長，而上式的二階估計模型可以表示為:

sj=arg min{E1(sj-1)+▽E1(sj-1)(s-sj-1)+

(21)

當(dāng)式(19)加入Tikhonov正則項，即:

E(s)=E1(s)+ρ‖▽s‖2+ρl‖s‖2

可以得到以下迭代公式:

sj=arg min{E1(sj-1)+▽E1(sj-1)(s-sj-1)+

(22)

忽略常數(shù)項，可以得到:

(23)

由于Tikhonov正則項是可分的，因此sj的計算可以變成解決每個元素的一維問題，它的迭代公式為:

sj=Tρtj(sj-1-tj▽E1(sj-1))

(24)

其中Tρ:Rn→Rn是收縮算子，定義為:

(25)

2.2 利用加速算子進行梯度估計

ISTA是梯度下降法的一種延伸，每次迭代只是利用當(dāng)前點的信息進行梯度估計，更新迭代點，算法迭代速度一般?；诖耍贗STA的基礎(chǔ)上引入加速算子，首先加速算子由當(dāng)前點和前一個點的線性組合構(gòu)成，然后利用加速算子重新進行梯度估計求出下一個迭代點，加速算子表示為：

yj+1=sj+aj(sj-sj-1)

(26)

其中:sj、sj-1分別代表當(dāng)前深度信息值和前一次深度信息值；sj-sj-1表示搜索方向；aj表示由當(dāng)前深度信息sj開始，沿著sj-sj-1所構(gòu)成的搜索方向進行迭代所需要的步長因子；而yj+1表示由當(dāng)前深度信息sj開始，沿著sj-sj-1所構(gòu)成的搜索方向進行步長為aj所得到的深度信息。利用加速算子重新進行梯度估計的示意圖如圖2所示。

圖2 加速算子梯度估計示意圖

在ISTA的基礎(chǔ)上，引入加速算子重新進行梯度估計求出下一個迭代點，因此式(24)可以寫成以下形式：

sj+1=Tρtj(yj+1-tj▽E1(yj+1))

(27)

2.3 割線線性搜索尋找最優(yōu)迭代步長

由于ISTA在迭代過程采用固定步長，優(yōu)化過程在靠近極小值點時，此時固定步長可能大于實際迭代步長，導(dǎo)致算法效率不佳；因此，采用線性搜索改變固定步長的限制。線性搜索主要有牛頓法、黃金分割法、割線法等，其中牛頓法收斂的條件是初始點充分靠近某一極值點，黃金分割法在求極值時需要事先確定搜索區(qū)間，而割線法能很快求出初始點附近的極值點，所以，采用割線線性搜索尋找最優(yōu)迭代步長，提高算法效率。

2.3.1 割線線性搜索

割線法的迭代公式表示如下：

(28)

其中:x(k)、x(k-1)是兩個初始點，x(k+1)是迭代后求出的點;f′(x(k-1))、f′(x(k))是x(k)、x(k-1)的導(dǎo)數(shù)值，其迭代過程如圖3所示。

圖3 割線線性搜索示意圖

2.3.2 尋找最優(yōu)迭代步長

為了求解最優(yōu)迭代步長，建立以下關(guān)于步長的一元函數(shù)：

t*=arg minψ(t)

(29)

其中：

ψ(t)=?(z(x,y,Δt)-I2(x,y))2dxdy+

ρ‖▽(s+td)‖2+ρl‖s+td‖2

(30)

其中：ψ(t)是關(guān)于t的一元函數(shù)，t是迭代步長，d是迭代方向，t*是最優(yōu)迭代步長，因此可以采用一維線性搜索方法來求極小值。

對于式(30)采用割線迭代格式表示為：

(31)

將利用割線線性搜索計算得到的最優(yōu)迭代步長t*代入式(27)，此時可以寫成以下形式：

(32)

FL-ISTA針對ISTA迭代步長固定的限制，引入割線線性搜索尋找最優(yōu)迭代步長，動態(tài)確定每次最優(yōu)迭代步長，該方法首先以大步長快速收斂到理想解附近，然后再以小步長逼近最優(yōu)解，從而提高算法效率。

3 算法流程

FL-ISTA步驟如下所示，算法流程如圖4所示。

步驟1 給出攝像機參數(shù)：焦距f,相距v,理想物距s0，凸透鏡半徑D,模糊度與模糊圓之間的常數(shù)λ，兩幅模糊圖像I1、I2，能量閾值ε，能量系數(shù)ρ，迭代步長t0，y1=s0,j=1，n1=α1=1, 0<αk<1。

步驟2 初始化深度信息為s0。

步驟3 計算式(12)，得到相對模糊度Δσ2。

步驟4 計算式(11)，得到擴散方程的解z(x,y,Δt)。

步驟5 用步驟4所得的解來計算能量式，如果能量小于閾值ε，則停止，此時的深度信息即為所求；否則，轉(zhuǎn)入步驟6。

步驟8 將步驟7所得到的yj+1代入式(27)。

步驟9 利用割線線性搜索求出最優(yōu)迭代步長t*，代入式(32)，求出此時的sj+1。

步驟10 將計算所得的sj+1代入式(14)，更新深度，回到步驟3，繼續(xù)迭代。

圖4 FL-ISTA流程

4 實驗結(jié)果及分析

為了驗證本文所提算法在求解離焦深度恢復(fù)動態(tài)優(yōu)化問題的正確性和有效性，分別采用ISTA、FISTA和MFISTA以及本文提出的FL-ISTA對標準500 nm尺度柵格進行深度信息重建實驗。標準納米尺度柵格高500 nm、寬1 500 nm。實驗中使用的是HIROX-7700顯微鏡，放大倍數(shù)為7 000，其余的攝像機參數(shù)為：焦距0.357 mm，理想物距3.4 mm，光圈大小為2。在配置為Intel i5-4900酷瑞雙核，主頻為3.3 GHz，內(nèi)存為8 GB的計算機上用Matlab 2012a對標準500 nm柵格進行實驗。

對標準500 nm尺度柵格的處理結(jié)果如圖5～8所示。其中：圖5是柵格的兩幅離焦圖像，圖5(a)是深度信息變化前的離焦圖像，圖5(b)是深度變化后的離焦圖像；圖6為標準500 nm尺度柵格的真實形貌；圖7是深度信息值收斂曲線；圖8(a)～(d)分別為ISTA、FISTA、MFISTA和FL-ISTA處理標準500 nm尺度柵格的3D形貌恢復(fù)結(jié)果。圖中，平面坐標的單位是像素，像素的大小為115.36 nm×115.36 nm，縱坐標的單位是mm。

圖5 柵格的離焦圖像

圖6 真實的柵格3D形貌

根據(jù)文獻[7]的實驗結(jié)果，將能量閾值ε設(shè)置為2，能量系數(shù)ρ設(shè)為0.2。在此條件下，ISTA、FISTA、MFISTA在迭代200次左右趨于收斂，而FL-ISTA在小于200次內(nèi)就已經(jīng)開始收斂，為了便于比較，設(shè)置迭代次數(shù)為250，如圖7為四種算法250次內(nèi)迭代深度信息值收斂曲線。由圖7可以看出，在迭代初期，四種算法都能以較快的速度收斂，屬于正?，F(xiàn)象。但是由于在離焦深度恢復(fù)動態(tài)優(yōu)化問題中，F(xiàn)L-ISTA在ISTA的基礎(chǔ)上，利用當(dāng)前點和前一個點的線性組合構(gòu)成加速算子，重新進行梯度估計，加快迭代速度；引入割線性搜索尋找最優(yōu)迭代步長，提高了算法效率。由60～120次的迭代過程中明顯可以看出，F(xiàn)L-ISTA得到的深度信息值下降最快，并且趨于收斂時得到的深度信息值相對較小，此時更加接近真實的標準500 nm柵格尺度，使得恢復(fù)出來的3D形貌尺度更精確。

圖7 四種算法深度信息值收斂曲線

由圖8(a)～(d)可以看出，四種算法重建的3D形貌，盡管在邊緣處誤差稍微大，但重建的3D形貌大致符合標準500 nm尺度柵格的整體形貌。從重建的3D形貌明顯看出，前三種算法重建的3D形貌中深度信息值相差不大，分別為600 nm、590 nm、590 nm左右，而FL-ISTA中深度信息值明顯下降很多，深度信息值為540 nm左右，更接近標準500 nm柵格尺度，使得恢復(fù)出來的3D形貌更精確。并且，在相同的條件下，四種算法對標準納米尺度柵格進行深度重建，經(jīng)多次實驗結(jié)果表明，ISTA、FISTA、MFISTA深度重建平均運行時間分別為105 s、105 s、100 s左右，而本文FL-ISTA的平均運行時間僅為30 s左右，一定程度上提高了算法的效率。

(33)

(34)

(35)

圖8 四種算法重建的3D形貌

由如圖9(a)～(d)可以看出，ISTA、FISTA、MFISTA所得的相對誤差曲面變化不大，最大相對誤差為80 nm、70 nm、70 nm左右，但是從FL-ISTA的相對誤差曲面來看，最大相對誤差只有20 nm左右，適合相對誤差較小的微納米操作場合。

圖9 四種算法的嘗試誤差曲面

由圖10可以看出，ISTA、FISTA和MFISTA均方差分別為0.05、0.048、0.045，而FL-ISTA均方差為0.041，與ISTA相比均方差下降了18個百分點，與FISTA和MFISTA相比均方差也較小，因此，F(xiàn)L-ISTA重建的微觀3D形貌尺度誤差相對穩(wěn)定。因為對于像細胞等樣品的觀測與測量，重建的3D形貌尺度與真實的3D形貌尺度誤差過大或過小，會導(dǎo)致對樣品尺度的估計不準確，進而可能造成操作過程中對樣品的破壞，從而造成嚴重后果。

由式(35)可以得出，ISTA、FISTA和MFISTA重建3D形貌平均誤差分別為161 nm、160 nm、158 nm，而FL-ISTA平均誤差僅為96 nm，與ISTA相比平均誤差下降了40個百分點，與FISTA和MFISTA相比平均誤差也較小，F(xiàn)L-ISTA重建3D形貌尺度與真實3D形貌尺度相差不大，可以利用重建3D形貌尺度來估計真實3D形貌尺度，以便于顯微鏡下對樣品的操作。

圖10 四種算法均方差的收斂曲線

最后為了更好地檢驗本文算法的性能，分別利用ISTA和FL-ISTA對于導(dǎo)電探針進行深度信息恢復(fù)實驗，如圖11為導(dǎo)電探針的離焦圖像，圖11(a)為深度變化前的離焦圖像，圖11(b)為深度變化后的離焦圖像；圖12為處理導(dǎo)電探針的實驗結(jié)果。

由圖12可以看出，F(xiàn)L-ISTA恢復(fù)的形貌與ISTA恢復(fù)的形貌有明顯的不同，并且在相同的條件下，經(jīng)過多次實驗可以得出，ISTA恢復(fù)形貌所需平均運行時間為180 s左右，而FL-ISTA恢復(fù)形貌所需平均運行時間僅為110 s，大大提升了算法的效率。

圖11 導(dǎo)電探針的離焦圖像

圖12 ISTA和FL-ISTA重建的導(dǎo)電探針3D形貌

5 結(jié)語

本文提出的FL-ISTA針對ISTA求解離焦深度恢復(fù)動態(tài)優(yōu)化問題迭代效率不佳，使得重建3D形貌精度不高的缺點進行優(yōu)化改進，利用加速算子重新進行梯度估計，更新迭代點；結(jié)合割線線性搜索尋找最優(yōu)迭代步長，提高算法效率?；跇藴?00 nm尺度柵格的實驗結(jié)果表明，與ISTA、FISTA、MFISTA相比，F(xiàn)L-ISTA的收斂速度更快，它恢復(fù)的3D形貌深度信息下降更快，相對誤差、均方差、平均誤差更小，更接近標準500 nm柵格尺度，有效地提高了求解離焦深度恢復(fù)動態(tài)優(yōu)化問題的效率和微觀3D形貌重建的精度。由于在迭代過程采用割線線性搜索來尋找最優(yōu)迭代步長仍需一段時間，所以接下來會進一步探索非線性搜索方法來提高算法效率。