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

李志勇，張 鵬

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

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

1 偏振圖像融合數學模型

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

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

2 融合快速算法分析及實現

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

2.1 總體優化

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

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

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

2.2 偏振度的快速運算

式(6)主要是進行平方根運算和除法運算，雖然可以采用標準算法，但運算周期較長，而且不利于硬件發揮并行運算優勢。

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

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

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

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

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

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

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

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

2.3 偏振角的快速運算

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

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

arctan(x)的Taylor級數展開式為式15。

當x＞0時， arctan(x)的第m項的截斷誤差為x2m+1/(2m+1)。如果m=1，就表示可以不經運算，直接用x作為arctan(x)的值。但在本應用中，還要滿足x2m+1/(2m+1)≤0.024 64(m=1)，從而可得，只有x≤0.419時，才能直接用x作為arctan(x)的值，簡化運算。

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

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

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

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

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

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

(3)最后，根據x的正負調整運算結果。當x≥0， θ=θ′；當x＜0，則有θ=π-θ′。

3 實驗測試及結論

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

實驗時，通過Camera Link總線控制三臺偏振相機同時曝光[12]，在得到偏振強度圖像后，觸發DSP中斷，然后DSP讀取圖像數據，進行配準等預處理后，開始融合處理。在融合處理中，針對本文提出的快速算法進行了DSP編程。限于篇幅，只列出了如圖2所示的典型過程。

圖1 基于DSP的偏振成像探測系統

圖2 偏振圖像融合快速運算典型過程

實驗中，分別對地物復雜、較復雜和簡單場景采集灰階偏振圖像，其像素分辨率為1024×1024。對同一次采集圖像， DSP系統分別用原始算法和本文算法進行融合處理，并通過DSP調試工具記錄兩種算法的處理時間，然后將融合圖像傳輸給計算機顯示，并進行數據分析。圖3 是對較復雜場景偏振圖像融合處理結果，其中從左到右依次為融合運算后得到的偏振強度圖、偏振度圖和偏振角圖。

圖3 偏振圖像融合處理結果

表1 復雜場景偏振圖像融合算法測試結果

表2 較復雜場景偏振圖像融合算法測試結果

表3 簡單場景偏振圖像融合算法測試結果

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

4 結束語

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

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

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

[ 3] 羅龍智，周南，羅海.整數開平方快速算法及其定點DSP實現[ J] .微計算機信息， 2007， 23(3-2)：157-158.

[ 4] 申伯純，周學軍.應用逐次遞進二分法在DSP上快速計算平方根[ 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] 伍微，劉小匯，李崢嶸， 等.實現定點DSP匯編層反正切函數的差分進化算法[ J].系統工程與電子技術， 2005， 27(5)：926-928.

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

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

[ 11] 楊偉鋒，洪津.無人機載偏振CCD相機光機系統設計[ J] .光學技術， 2008， 34(3)：469-472.

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