基于DSP的偏振圖像快速融合研究*

李志勇，張 鵬

(解放軍炮兵學(xué)院， 合肥230031)

由于能夠提供比普通的強(qiáng)度圖像和光譜圖像更加豐富的目標(biāo)場(chǎng)景信息，可見光偏振探測(cè)在相關(guān)社會(huì)和軍事領(lǐng)域日益受到重視。偏振探測(cè)的關(guān)鍵是獲取能反映研究對(duì)象特征的光偏振態(tài)信息，但必須先采集多個(gè)起偏角的偏振圖像，再通過像素級(jí)融合運(yùn)算才能得到[1]。雖然該運(yùn)算不復(fù)雜，但為了保證良好的識(shí)別分析能力，偏振圖像的分辨率一般要求較高，像素多，這使得總運(yùn)算量很大，耗時(shí)長(zhǎng)，在實(shí)際應(yīng)用中受限。目前國(guó)內(nèi)外對(duì)可見光偏振探測(cè)的研究主要集中在圖像特征分析上，而對(duì)圖像的快速融合在包括參考文獻(xiàn)[2]的公開文獻(xiàn)上基本沒有提及。因此，本文從快速化的角度改進(jìn)該運(yùn)算算法，以提高偏振探測(cè)系統(tǒng)的實(shí)時(shí)性。

1 偏振圖像融合數(shù)學(xué)模型

在偏振探測(cè)中，常用斯托克斯(Stokes)矢量表示法來描述光偏振態(tài)，包括I(光波總強(qiáng)度)、Q(水平方向上線偏振光強(qiáng)度)、U(45°/135°方向上線偏振光強(qiáng)度)和V(圓偏振光強(qiáng)度)。在實(shí)際應(yīng)用中， V往往非常小，可以忽略不計(jì)。為了得到I、Q和U，需要獲取不同起偏角的偏振圖像， 然后通過式(1)～式(5)的數(shù)學(xué)模型[1]，進(jìn)行像素級(jí)的融合運(yùn)算，得到能反映目標(biāo)特性的線偏振度P和偏振角θ。

其中， I0°、I60°和I120°分別表示起偏角0°， 60°和120°的偏振強(qiáng)度值。此外，式(4)和式(5)還有約束條件，詳見參考文獻(xiàn)[1]。

2 融合快速算法分析及實(shí)現(xiàn)

通過對(duì)其它快速算法及其工程實(shí)踐分析可知，要提高算法的實(shí)時(shí)性，根本上來說，就是在保證對(duì)運(yùn)算結(jié)果影響很小的前提下，盡量增加算法中簡(jiǎn)單而容易實(shí)現(xiàn)的運(yùn)算的比例，如加法和乘法。同時(shí)還要考慮利用DSP處理系統(tǒng)的運(yùn)算特點(diǎn)，如定點(diǎn)運(yùn)算速度快于浮點(diǎn)運(yùn)算。

2.1 總體優(yōu)化

分析偏振圖像融合數(shù)學(xué)模型可知，如果按其原步驟依次運(yùn)算，基本的浮點(diǎn)運(yùn)算較多，包括6次乘法、1次加法和兩次除法。分析可知，造成此現(xiàn)象的根本原因是I、Q和U都是浮點(diǎn)數(shù)，使得后續(xù)運(yùn)算都只能是浮點(diǎn)型。從優(yōu)先考慮整數(shù)運(yùn)算的角度出發(fā)，設(shè)I′=I0°+I60°+I120°， Q′=2I0°-I60°-I120°， U′=I60°-I120°。然后將式(1)、式(2)和式(3)代入式(4)和式(5)，則有式(6)和式(7)。并可再設(shè)整數(shù)SO=Q′2+3U′2

對(duì)原步驟優(yōu)化后，按照先算整數(shù)I′、Q′和U′，再算SO、式(6)和式(7)的過程，浮點(diǎn)運(yùn)算減少了5次。而且加法和乘法的總次數(shù)也減少了1次。

在進(jìn)行上述優(yōu)化后，式(6)和式(7)，即偏振度和偏振角運(yùn)算成為最復(fù)雜的部分。需要對(duì)它們進(jìn)一步分析。

2.2 偏振度的快速運(yùn)算

式(6)主要是進(jìn)行平方根運(yùn)算和除法運(yùn)算，雖然可以采用標(biāo)準(zhǔn)算法，但運(yùn)算周期較長(zhǎng)，而且不利于硬件發(fā)揮并行運(yùn)算優(yōu)勢(shì)。

不少文獻(xiàn)，如文獻(xiàn)[3 -4]，提出了一些結(jié)合實(shí)踐的快速開平方根算法，其中成熟且較合適的是以加法和乘法為主的牛頓迭代法。其公式為：

式(8)中SO表示需要開平方根的數(shù)， Sn表示第n次得到的平方根值。該方法的首要關(guān)鍵是估計(jì)合適的第一個(gè)近似根S1。如果S1與真實(shí)根越接近，那么需要的迭代次數(shù)就越少。對(duì)此，可以借鑒參考文獻(xiàn)[5]，引入Taylor級(jí)數(shù)法，并利用計(jì)算機(jī)中浮點(diǎn)數(shù)的定點(diǎn)表示方法，將S1的浮點(diǎn)數(shù)估計(jì)轉(zhuǎn)換為容易確定的32位整數(shù)估計(jì)問題，并使式(8)只含加法和乘法。但通過實(shí)踐發(fā)現(xiàn)， Sn的倒數(shù)運(yùn)算過于簡(jiǎn)化，會(huì)影響最終結(jié)果精度。因此，本文對(duì)此進(jìn)行了改進(jìn)，引入了牛頓迭代法求倒數(shù)公式，初步得到如下算法過程。

(1)設(shè)整數(shù)CI0和單精度浮點(diǎn)數(shù)S1共用同一存儲(chǔ)地址，整數(shù)CI1和單精度浮點(diǎn)數(shù)S′1(S1的倒數(shù))共用同一存儲(chǔ)地址，按式(9)初始運(yùn)算。其中，因?yàn)镮′是整數(shù)，其取值范圍與像素精度有關(guān)，有限且預(yù)知，可以事先算出它的平方的倒數(shù)，然后在實(shí)時(shí)處理中用查表法得到。

(2)由于分別共用同一地址，可用單精度浮點(diǎn)數(shù)S1和S′1中的十六進(jìn)制數(shù)值表示整數(shù)CI0和CI1，并分別用式(10)和式(11)分別變換CI0和CI1，以得到合適的S1[2]：

(3)整數(shù)CI0和CI1中的數(shù)值轉(zhuǎn)而表示單精度浮點(diǎn)數(shù)S1和S′1。此時(shí)，用式12(牛頓迭代法求倒數(shù)公式)得到S′1：

(4)將S1和S′1代入式13，進(jìn)行第一次迭代求出平方根：

(5)檢驗(yàn)兩次平方根的差是否滿足精度要求。如果不滿足，繼續(xù)下一次迭代運(yùn)算。

此時(shí)，需要解決的是如何加快收斂速度問題。對(duì)此，文獻(xiàn)[6]等提出了一些具有三階收斂速度的方法。而在本應(yīng)用中，有限的像素精度(原始圖像為12 bit，融合圖像為8 bit)決定了運(yùn)算精度是有限的，可以據(jù)此減少迭代次數(shù)，加快收斂速度。在物理上， P表示的是從非偏振光情況下的0到全偏振光情況下的1之間的浮點(diǎn)數(shù)。在融合圖像中，每個(gè)P都要用一個(gè)8 bit像素表示，也即P要轉(zhuǎn)換為區(qū)間[0， 255]中的偏振度圖像素值PL，即有PL=255×P。設(shè)平方根的運(yùn)算誤差為浮點(diǎn)數(shù)ES， PL的運(yùn)算誤差為整數(shù)EP，則根據(jù)式(6)有EP=255×ES/I′。在保證EP不大于1個(gè)像素灰度值的要求下，必須有ES≤I′/255，則ES的最大允許誤差值為I′/255。根據(jù)式1、2和3可知， SO與I′的大小是相關(guān)的。在牛頓迭代法中，雖然ES隨著SO增大而增加，但同時(shí)I′也在增加，仍然有可能保證ES≤I′/255。在本應(yīng)用中， SO與I′的取值范圍是有限且已知的。因此，通過軟件進(jìn)行了測(cè)試，代入SO與I′的所有取值，得到牛頓迭代法求平方根的運(yùn)算結(jié)果，并與C語言庫(kù)的平方根函數(shù)的運(yùn)算結(jié)果比較，得到ES。結(jié)果表明只需要一次迭代，就能滿足精度要求。

2.3 偏振角的快速運(yùn)算

與式(9)類似，式(7)中的Q′和U′的倒數(shù)也能提前算出，采用查表法變除法為乘法。其關(guān)鍵還是要解決反正切的快速運(yùn)算問題。不少文獻(xiàn)都提供了一些反正切快速算法，其中，文獻(xiàn)[ 7]和[8]的算法還是較復(fù)雜，而文獻(xiàn)[9]雖然提出了一種較好的查找表法，但沒有利用運(yùn)算精度有限這個(gè)條件。

設(shè)arctan(x)的誤差為EA， 則T的誤差 ET為40.584 5×EA。為了不影響融合圖像精度，需要保證ET≤1(最小像素灰度值)，從而有EA≤0.024 64。

arctan(x)的Taylor級(jí)數(shù)展開式為式15。

當(dāng)x＞0時(shí)， arctan(x)的第m項(xiàng)的截?cái)嗾`差為x2m+1/(2m+1)。如果m=1，就表示可以不經(jīng)運(yùn)算，直接用x作為arctan(x)的值。但在本應(yīng)用中，還要滿足x2m+1/(2m+1)≤0.024 64(m=1)，從而可得，只有x≤0.419時(shí)，才能直接用x作為arctan(x)的值，簡(jiǎn)化運(yùn)算。

如果x＞0.419，經(jīng)分析可知， x越大，需要展開的項(xiàng)數(shù)也越多，運(yùn)算量增加，不宜再用Taylor級(jí)數(shù)法。此時(shí)可以利用反正切函數(shù)的周期性和對(duì)稱性進(jìn)行運(yùn)算。其思路是將x轉(zhuǎn)換為取值區(qū)間是[ 0， 1]的x′，從而求出arctan(x)在對(duì)應(yīng)區(qū)間[ 0， π]中的值。其過程如下：

(1)分別求取U′和Q′的絕對(duì)值Uf′和Qf′，保證它們都是正數(shù)，便于后續(xù)運(yùn)算；

(2)如果Uf′≤Qf′，則有x′=Uf′/Qf′， θ′=arctanx′；否則， x′=Qf′/)， θ′=π/2-arctanx′。

根據(jù)式(11)， θ的區(qū)間[0， π]被255等分，每一等分弧度α=π/255。當(dāng)有kα≤θ＜(k+1)α (k=0， 1，2， 3…， 255)，使得T=k。這表明，在[kα，(k+1)α]區(qū)間內(nèi)的任何θ所得的T都是相同的。進(jìn)一步可得：tan(kα)≤tanθ＜tan((k+1)α)。因?yàn)檎泻瘮?shù)和反正切函數(shù)都是正比例函數(shù)，所以，與θ運(yùn)算特點(diǎn)一致，在區(qū)間[tan(kα)， tan((k+1)α)]內(nèi)的任意值運(yùn)算所得的T都等于通過tan(kα)運(yùn)算所得的T。這表明，在本應(yīng)用條件下，沒必要求出所有值，只要事先求出有限點(diǎn)的反正切值，通過判斷需要進(jìn)行反正切運(yùn)算的值的區(qū)間范圍[tan(kα)， tan((k+1)α)]，就能用查表法直接賦值。同理，對(duì)于arctanx′，每一等分弧度α′=2α，根據(jù)arctanx′∈[0， π/4]，可推導(dǎo)出 x′的數(shù)值范圍判斷區(qū)間為[tan(kα′)， tan((k+1)α′)] (k=0， 1...31)。因?yàn)樾枰蟮膖ankα′≥0.419， 從而得k≥16.1。所以，從tan16α′開始，查找表中需要建立的數(shù)值點(diǎn)只有16個(gè)。

使用查找表法，存在如何快速查找的問題。如果先判斷再查找，采用折半法搜索，最多需要4 次。搜索次數(shù)不多，但在加入多次判斷的情況下，由于跳轉(zhuǎn)較多，不利于并行運(yùn)算。分析正切值與各段區(qū)間的關(guān)系及規(guī)律發(fā)現(xiàn)，當(dāng)有式(9)時(shí)， Bn(整數(shù)n的取值范圍為[16， 31])的整數(shù)部分基本上形成一個(gè)差值為1的等差數(shù)列，其范圍為0-20，正好可以對(duì)應(yīng)作為數(shù)組下標(biāo)。數(shù)組的內(nèi)容為kα′。

但是，通過計(jì)算可知， Bn的取值范圍是不連續(xù)的，缺少7、12、16和18等值。需要把這些值代入式12中，反推得到對(duì)應(yīng)的反正切值，并賦予相應(yīng)數(shù)組內(nèi)容。這樣，在式(12)中，用x′代替tankα′進(jìn)行運(yùn)算，用一次乘法和兩次加法，就能得到x′對(duì)應(yīng)的數(shù)組下標(biāo)，進(jìn)而查表取值。相比文獻(xiàn)[6]中需要建立256 byte的數(shù)值表，本文方法只需要42 byte，而且查找速度更快。

(3)最后，根據(jù)x的正負(fù)調(diào)整運(yùn)算結(jié)果。當(dāng)x≥0， θ=θ′；當(dāng)x＜0，則有θ=π-θ′。

3 實(shí)驗(yàn)測(cè)試及結(jié)論

在以DSP系統(tǒng)為核心的偏振成像探測(cè)系統(tǒng)樣機(jī)[10-11]上進(jìn)行了測(cè)試。如圖1所示。各臺(tái)偏振相機(jī)由圖像傳感器與偏振片和延遲器組成，分別只能接收某一起偏角的偏振強(qiáng)度圖像。 DSP系統(tǒng)選擇了TI公司的TMS320C6711D為核心芯片，其工作頻率最高為250 MHz。計(jì)算機(jī)用于實(shí)驗(yàn)控制及數(shù)據(jù)分析。

實(shí)驗(yàn)時(shí)，通過Camera Link總線控制三臺(tái)偏振相機(jī)同時(shí)曝光[12]，在得到偏振強(qiáng)度圖像后，觸發(fā)DSP中斷，然后DSP讀取圖像數(shù)據(jù)，進(jìn)行配準(zhǔn)等預(yù)處理后，開始融合處理。在融合處理中，針對(duì)本文提出的快速算法進(jìn)行了DSP編程。限于篇幅，只列出了如圖2所示的典型過程。

圖1 基于DSP的偏振成像探測(cè)系統(tǒng)

圖2 偏振圖像融合快速運(yùn)算典型過程

實(shí)驗(yàn)中，分別對(duì)地物復(fù)雜、較復(fù)雜和簡(jiǎn)單場(chǎng)景采集灰階偏振圖像，其像素分辨率為1024×1024。對(duì)同一次采集圖像， DSP系統(tǒng)分別用原始算法和本文算法進(jìn)行融合處理，并通過DSP調(diào)試工具記錄兩種算法的處理時(shí)間，然后將融合圖像傳輸給計(jì)算機(jī)顯示，并進(jìn)行數(shù)據(jù)分析。圖3 是對(duì)較復(fù)雜場(chǎng)景偏振圖像融合處理結(jié)果，其中從左到右依次為融合運(yùn)算后得到的偏振強(qiáng)度圖、偏振度圖和偏振角圖。

圖3 偏振圖像融合處理結(jié)果

表1 復(fù)雜場(chǎng)景偏振圖像融合算法測(cè)試結(jié)果

表2 較復(fù)雜場(chǎng)景偏振圖像融合算法測(cè)試結(jié)果

表3 簡(jiǎn)單場(chǎng)景偏振圖像融合算法測(cè)試結(jié)果

在計(jì)算機(jī)中，以原始算法的融合圖像為標(biāo)準(zhǔn)，用均方差MSE來評(píng)價(jià)本文算法所得融合圖像的質(zhì)量。最后，分別得到三種場(chǎng)景的50次測(cè)試比較結(jié)果的平均值如表1、表2和表3所示，與前者相比，快速算法明顯地提高了處理實(shí)時(shí)性，而對(duì)圖像質(zhì)量影響很小。同時(shí)可知，由于偏振圖像融合處理只與不同偏振方向強(qiáng)度差別有關(guān)，因此不同場(chǎng)景復(fù)雜程度對(duì)處理時(shí)間影響很小。此外需要說明的是，表1中處理時(shí)間包括圖像數(shù)據(jù)存取時(shí)間(兩種算法相同)。

4 結(jié)束語

為了提高偏振圖像融合處理的速度，本文通過理論算法和硬件資源結(jié)合進(jìn)行了快速算法研究，并在基于DSP的實(shí)驗(yàn)平臺(tái)測(cè)試，取得了較好的效果。由于所研制的偏振探測(cè)系統(tǒng)只是樣機(jī)，融合處理速度還有較大的提高余地。例如下一步擬采用硬件性能更好的單個(gè)或多個(gè)DSP芯片，以及增強(qiáng)程序優(yōu)化程度。

[ 1] 王新，王學(xué)勤，孫金祚.基于偏振成像和圖像融合的目標(biāo)識(shí)別技術(shù)[ J] .激光與紅外， 2007， 37(7)：676-678.

[ 2] 張晶晶，方勇華，陳曉寧.基于小波的偏振圖像融合算法及性能評(píng)價(jià)[ J] .合肥工業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版)， 2009， 32(7)：1101-1105.

[ 3] 羅龍智，周南，羅海.整數(shù)開平方快速算法及其定點(diǎn)DSP實(shí)現(xiàn)[ J] .微計(jì)算機(jī)信息， 2007， 23(3-2)：157-158.

[ 4] 申伯純，周學(xué)軍.應(yīng)用逐次遞進(jìn)二分法在DSP上快速計(jì)算平方根[ J] .電子工程師， 2008， 33(4)：59-61.

[ 5] Chris Lomont.Fast Inverse Square Root[ EB/OL] .www.lomont.org/math/2003.2003.

[ 6] OZban A Y.Some Variants of Newton' s Methods[ J].Applied Mathematics Letters， 2004， 17：677-682.

[ 7] Dave Van Ess.Algorithm-ArcTan as Fast as You Can[ EB/OL].www.cypress.com/? rID=2701.2006.

[ 8] 伍微，劉小匯，李崢嶸， 等.實(shí)現(xiàn)定點(diǎn)DSP匯編層反正切函數(shù)的差分進(jìn)化算法[ J].系統(tǒng)工程與電子技術(shù)， 2005， 27(5)：926-928.

[ 9] 劉禮剛，曾延安，常大定.基于FPGA的反正切函數(shù)的優(yōu)化算法[ J] .微計(jì)算機(jī)信息， 2007， 23(6-2)：203-204.

[ 10] 王峰，洪津，喬延利，等.專用偏振成像智能遙感器設(shè)計(jì)研究[ J] .傳感技術(shù)學(xué)報(bào)， 2007， 20(6)：1448-1452.

[ 11] 楊偉鋒，洪津.無人機(jī)載偏振CCD相機(jī)光機(jī)系統(tǒng)設(shè)計(jì)[ J] .光學(xué)技術(shù)， 2008， 34(3)：469-472.

[ 12] 李志勇， 袁魏華， 楊振華.基于TMS320C6711的Camera Link相機(jī)控制的實(shí)現(xiàn)[ J] .電子器件， 2007， 30(3)：972-975.