基于二次成像的光學(xué)圖像相位恢復(fù)

婁聯(lián)堂，但 威，陳佳騏，吳高林

(中南民族大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)學(xué)院，武漢 430074)

相位恢復(fù)(PR)技術(shù)是利用已知的強(qiáng)度信息，通過迭代找到最符合真實(shí)場強(qiáng)數(shù)據(jù)的相位分布.相位恢復(fù)算法在信號(hào)復(fù)原、設(shè)計(jì)光學(xué)衍射元件等方面有所應(yīng)用.迄今為止，已有很多有影響的相位恢復(fù)算法. 20世紀(jì)70年代Gerchbery和Saxton兩位學(xué)者提出了有實(shí)用意義的迭代相位恢復(fù)算法—GS算法[1]，即由已知像平面和衍射平面上的強(qiáng)度分布迭代算出相位值.1973年Misell仿照仿照GS算法方案，提出如何由已知兩張具有不同離焦量的離焦像來導(dǎo)出波函數(shù)的相位值.Fienup于1982年提出了誤差減小算法(Error Reduction algorithm, ER)[2]和混合輸入輸出算法(Hybrid Input-Output, HIO)，使得收斂速度大大提高.國內(nèi)學(xué)者楊國楨和顧本源在20世紀(jì)90年代提出了一個(gè)有效的算法，即楊-顧算法[3-6]，這個(gè)算法是根據(jù)輸入和輸出兩個(gè)面上的強(qiáng)度信息計(jì)算出輸入平面上的相位信息.

本文在楊-顧算法的基礎(chǔ)上，提出了一種基于二次成像光學(xué)圖像相位恢復(fù)的新算法.首先通過建立簡化的光學(xué)成像系統(tǒng)，得到光學(xué)圖像之間的薛定諤變換[7]關(guān)系，然后利用圖像薛定諤變換和楊-顧算法來恢復(fù)相位.

1 楊-顧算法

楊-顧算法原理[3]簡述如下：成像系統(tǒng)中任意的某兩平面上的波函數(shù)：

U1(x1)=|U1(x1)|ejφ1(x1),

U2(x2)=|U2(x2)|ejφ2(x2),

是通過幺正算符H(H+H=I，幺正算符即H+=H-1)聯(lián)系起來，即：

U2(x2)=HU1(x1)或|U2(x2)|=

e-jφ2(x2)HU1(x1),

U1(x1)=H-1U2(x2)或

|U1(x1)|=e-jφ1(x1)H-1U2(x2).

楊-顧算法恢復(fù)相位，就是要找到這樣的φ1和φ2，使得e-jφ2H·[|U1|ejφ1]盡可能的逼近|U2|.用距離D表示二者的逼近程度：

D2=‖|U2|-e-jφ2HU1‖2=

2ReTr[e-jφ1|U1|+HU2]}.

其中|U1|，|U2|表示N維列向量，e-jφ1和e-jφ2表示相位元素組成的N×N對(duì)角矩陣，Tr表示矩陣的跡，Re為復(fù)數(shù)的實(shí)數(shù)部分.

利用矩陣跡的輪換不變性，相位分布可以通過解下面的變分方程得到，即：

δφ1φ2‖|U2|-e-jφ2HU1‖=0.

解得：

ejφ2=Hejφ1|U1|/abs(Hejφ1|U1|),

(1)

ejφ1=H+ejφ2|U2|/abs(H+ejφ2|U2|),

(2)

已知振幅|U1|和|U2|，可以由(1)、(2)式迭代算出相位φ1和φ2.

2 光學(xué)成像系統(tǒng)模型

成像過程的建模是計(jì)算機(jī)視覺應(yīng)用的一個(gè)重要論題，如移動(dòng)的結(jié)構(gòu)、物體識(shí)別和位姿評(píng)估等等.在幾何學(xué)上，一張照片數(shù)據(jù)將三維實(shí)體空間變成二維圖像空間.用于計(jì)算機(jī)視覺的相機(jī)模型一般是立體投影模型，為了重建光學(xué)圖像的強(qiáng)度測量過程，可以采用一個(gè)簡化相機(jī)模型來實(shí)現(xiàn)，如圖1所示(與文獻(xiàn)[8]中的Figure 10-28類似).圖1給出了一個(gè)簡單的相機(jī)模型，相機(jī)的光圈置于z=z3平面，光源結(jié)合面(z=z5)位于凸透鏡鏡片的右邊，一般來說，相機(jī)成像的CCD面(z=z4或z=z6)會(huì)在結(jié)合面的左邊或右邊. 用f表示有效焦距 (或簡單焦距).用dl表示半徑，f滿足如下關(guān)系：

(3)

圖1 光學(xué)成像系統(tǒng)Fig.1 Optical imaging system

由文獻(xiàn)[8]可知，相機(jī)成像過程中，物體u1(x,y)經(jīng)過鏡片、光圈等的一系列變換處理聚集到CCD上形成密度圖像：

u4(x,y)=u1(x,y)**h12(x,y)**

tl(x,y)**h23(x,y)**t3(x,y)**h34(x,y)，

其中“**”表示兩個(gè)二元函數(shù)間的二維卷積.對(duì)上式作傅里葉變換有：

U4(ξ,η)=U1(ξ,η)H12(ξ,η)Tl(ξ,η)H23(ξ,η)T3(ξ,η)H34(ξ,η),

其中：

H23(ξ,η)=λB23z23jq(ξ,η;-λz23),

H34(ξ,η)=λB34z34jq(ξ,η;-λz34),

那么就有：

其中q(x,y;a)是二次相位信號(hào)函數(shù)，即q(x,y;a)=eja(x2+y2)，q*(x,y;a)為其共軛.q(x,y;a)的傅里葉變換函數(shù)是：

(4)

(5)

(6)

(6)式涵蓋了物距z12、像距z23+z34及焦距f三種變化情況，基本覆蓋了大多數(shù)光學(xué)成像情況.對(duì)(6)式作傅里葉逆變換得：

cF-1(ejt(ξ2+η2)F(u4(x,y))),

(7)

其中t為薛定諤變換尺度參數(shù)[9].

3 相位恢復(fù)算法

3.1 模型的建立

設(shè)用同一個(gè)相機(jī)在手動(dòng)調(diào)焦模式下采用不同的焦距獲取的同一目標(biāo)的兩幅圖像為|U1(x,y)|和|U2(x,y)|，現(xiàn)在可以假設(shè)連續(xù)信號(hào)分別是空間受限和帶寬受限的.圖像|U1(x,y)|和|U2(x,y)|的取樣點(diǎn)數(shù)目均為N=mn，m，n表示圖像|U1(x,y)|和|U2(x,y)|的寬和高.利用公式(7)并考慮到數(shù)字圖像的歸一化，那么圖像|U1(x,y)|和|U2(x,y)|滿足以下關(guān)系：

(8)

其中c,d是線性歸一化常數(shù)(同時(shí)考慮到(7)式中的常數(shù)c)，其作用是修正兩次成像得到的圖像中心及圖像大小的偏移，Ui(x,y)=|Ui(x,y)|ejφi,i=1,2.

根據(jù)(8)式以及|U1(x,y)|和|U2(x,y)|來恢復(fù)丟失的相位信息φ1(x,y)和φ2(x,y)，|U2(x,y)|與c·e-jφ2S(U1(x,y))+d·1的逼近程度，引入“距離”D2來描述，D2的定義[5,6]為：

D2=‖c·|U2|+d·1-e-jφ2F-1e-jt(ξ2+η2)Fejφ1|U1|‖2=‖U3-e-jφ2HU1‖2,

(9)

其中|U1|，|U2|和1都是N維列向量，|U1|，|U2|由圖像|U1(x,y)|和|U2(x,y)|的離散數(shù)值矩陣?yán)背蒒維列向量.e-jΦ2為N×N對(duì)角矩陣，U3=c·|U2|+d·1,F-1,F分別為圖像的傅里葉逆變換和正變換的變換矩陣，圖像薛定諤變換的變換矩陣H=F-1e-jt(ξ2+η2)F為幺正算符.

3.2 模型求解及恢復(fù)算法

由(9)式可得：

(10)

那么相位恢復(fù)問題可以轉(zhuǎn)化為求解下列變分方程：

δφ1(D2)=0,δφ2(D2)=0.

而圖像歸一化的系數(shù)c,d則可以由式求導(dǎo)算出.

于是，

ejφ2=Hejφ1|U1|abs(Hejφ1|U1|).

(11)

ejφ1=H+ejφ2U3/abs(H+ejφ2U3).

(12)

c={NRe[Tr(e-jφ2HU1|U2|+)]-

Tr(|U2|1+)Re[Tr(e-jφ2HU11+)]}/

{NTr(|U2||U2|+)-(Tr(|U2|1+))2},

(13)

d={Tr(|U2||U2|+) Re[Tr(e-jφ2HU11+)]-

Tr(|U2|1+)Re[Tr(e-jφ2HU1|U2|+)]}/

{NTr(|U2||U2|+)-(Tr(|U2|1+))2},

(14)

其中abs表示取絕對(duì)值，e-jφ1為對(duì)角矩陣.

下面給出基于二次成像的光學(xué)圖像相位恢復(fù)的算法步驟：

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

為了驗(yàn)證新的相位恢復(fù)算法，考慮了室內(nèi)和室外兩種情況，相機(jī)選用Nikon DX40單反數(shù)碼相機(jī)，實(shí)驗(yàn)是在Windows7下使用MATLAB7.10.0下完成的.

采用Nikon DX40單反數(shù)碼相機(jī)(手動(dòng)調(diào)焦模式)，通過調(diào)節(jié)相機(jī)焦距對(duì)室內(nèi)水杯進(jìn)行拍照取樣，其ISO曝光度均為100.本文選取了焦距分別為18mm和25mm兩張光學(xué)圖像(圖2)，經(jīng)縮小后的圖像大小均為581×389.因此，在上述迭代算法中，N的值就是581×389，圖2(a)為|U1(x,y)|，圖2(b)為|U2(x,y)|，|U1|，|U2|是由581×389矩陣?yán)背蒒=581×389維列向量，而ejφ1，ejφ2為N×N對(duì)角矩陣.在實(shí)驗(yàn)中，利用圖像的薛定諤變換和MATLAB的矩陣點(diǎn)乘來實(shí)現(xiàn)相位恢復(fù)的快速算法.初始相位φ1，φ2均取為581×389的零矩陣.薛定諤變換尺度參數(shù)t設(shè)為0.05，誤差限ESP設(shè)為100(即平均每個(gè)像素的誤差小于100/N).圖3給出了通過恢復(fù)的相位重建第二幅圖像的過程，隨著迭代次數(shù)的增加，重建誤差越來越小，其中重建誤差是根據(jù)(9)式計(jì)算得到，圖3(a)、(b)、(c)的誤差ERR分別為49600、9600、100.

(a) f=18 mm (b) f=25 mm圖2 室內(nèi)場景同一目標(biāo)不同焦距下二次成像Fig.2 Two imaging picture of indoor scene with different focal length

(a) ERR=49600 (b) ERR=9600 (c) ERR=100圖3 通過相位恢復(fù)重建的圖像及重建誤差Fig.3 Reconstructed images and errors by phase retrieval

圖4(a)、(b)是用同一臺(tái)相機(jī)對(duì)室外水杯分別選取焦距為18mm和25mm的二次成像結(jié)果，圖4(c)是利用圖4(a)的灰度、恢復(fù)的相位數(shù)據(jù)及(8)式重建的圖像，比較圖4(b)、(c)可以看出重建的效果和精度還是比較高的.

(a) f=18mm (b) f=25mm (c) 通過相位重建后的圖像圖4 室外場景二次成像及通過相位重建后的圖像Fig.4 Two imaging pictures of outdoor scene and reconstructed images by phase retrieval

圖5給出了圖2(a)、(b)及圖4(a)、(b)最終相位恢復(fù)結(jié)果.其中相位圖是直接將0到2π范圍內(nèi)的相位角歸一化到0～255來實(shí)現(xiàn)的，因此會(huì)出現(xiàn)一定的誤差.從相位圖中能大致看出杯子、凳子及背景分界線的形狀，因此表明恢復(fù)的相位圖有一定的可信度.

(a)圖2(a)相位圖 (b)圖2(b)相位圖 (c)圖4(a)相位圖 (d)圖4(b)相位圖 圖5 相位恢復(fù)結(jié)果Fig.5 Results of phase retrieval

實(shí)驗(yàn)中，為了不分散對(duì)算法的注意力，文中假設(shè)薛定諤變換尺度參數(shù)t為某一固定值(實(shí)驗(yàn)中取為0.05)，從恢復(fù)圖像以及得到的相位圖像可以看出恢復(fù)的效果比較好.但在實(shí)際中，這個(gè)參數(shù)是不確定的，它與相機(jī)成像模型參數(shù)有關(guān)，如何求出參數(shù)t，將是下一步工作的重點(diǎn)之一.