基于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.