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

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

0 引言

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

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

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

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

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

1.1 離焦成像模型

在圖像處理中，圖像的模糊化通?？梢杂镁矸e的形式表示如下：

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

(1)

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

(2)

并且:

σ2=2tc

(3)

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

(4)

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

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

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

(5)

(6)

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

(7)

并且在Δt時(shí)刻，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分別為攝像機(jī)的像距、焦距和物距，D為凸透鏡半徑。

1.2 深度信息重建

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

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

(11)

其中相對(duì)模糊度又可以表示為:

(12)

圖1 深度信息原理圖

本文定義:

(13)

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

(14)

為了更好地解決離焦深度恢復(fù)動(dòng)態(tài)優(yōu)化問(wèn)題，建立如下目標(biāo)函數(shù)來(lái)計(jì)算熱輻射方程:

(15)

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

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

(16)

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

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

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

(17)

因此，可以通過(guò)求解如下優(yōu)化問(wèn)題來(lái)求解景物深度信息:

(18)

s.t. (11)和(14)

2 FL-ISTA

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

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

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

(19)

于是，對(duì)于無(wú)約束最優(yōu)化問(wèn)題:

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

解決上式的梯度算法為:

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

(20)

其中tj表示迭代步長(zhǎng)，而上式的二階估計(jì)模型可以表示為:

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

(21)

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

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

可以得到以下迭代公式:

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

(22)

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

(23)

由于Tikhonov正則項(xiàng)是可分的，因此sj的計(jì)算可以變成解決每個(gè)元素的一維問(wèn)題，它的迭代公式為:

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

(24)

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

(25)

2.2 利用加速算子進(jìn)行梯度估計(jì)

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

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

(26)

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

圖2 加速算子梯度估計(jì)示意圖

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

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

(27)

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

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

2.3.1 割線線性搜索

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

(28)

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

圖3 割線線性搜索示意圖

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

為了求解最優(yōu)迭代步長(zhǎng)，建立以下關(guān)于步長(zhǎng)的一元函數(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是迭代步長(zhǎng)，d是迭代方向，t*是最優(yōu)迭代步長(zhǎng)，因此可以采用一維線性搜索方法來(lái)求極小值。

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

(31)

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

(32)

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

3 算法流程

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

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

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

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

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

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

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

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

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

圖4 FL-ISTA流程

4 實(shí)驗(yàn)結(jié)果及分析

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

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

圖5 柵格的離焦圖像

圖6 真實(shí)的柵格3D形貌

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

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

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

(33)

(34)

(35)

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

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

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

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

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

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

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

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

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

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

5 結(jié)語(yǔ)

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