高動態(tài)范圍視頻的多曝光圖像序列快速融合

樸永杰，徐 偉，王紹舉，陶淑蘋

（中國科學院 長春光學精密機械與物理研究所 小衛(wèi)星技術國家地方聯(lián)合工程研究中心，吉林 長春 130033）

1 引 言

普通相機在拍攝大動態(tài)范圍場景時很難在一次曝光圖像中展現(xiàn)所有的細節(jié)。因此人們開發(fā)了各種基于曝光時間控制的自動曝光方法［1－3］，以捕獲更多的場景細節(jié)。然而這些方法只能兼顧圖像的一個方面：畫面亮度高的部分可以正常曝光而亮度低的部分不能正常曝光，或者亮度低的部分正常曝光而亮度高的部分曝光過度。為此，文獻［4］首次從多幅LDR圖像（8～10幅）推算出系統(tǒng)的輻照圖，進而將多幅LDR圖像合成為一幅高動態(tài)范圍（HDR）圖像。

此后出現(xiàn)了眾多多次曝光法合成HDR圖像的研究［5－8］，但大多針對單幅圖像的 HDR 合成，沒有考慮HDR視頻獲取的實時性問題。文獻［9］對多曝光圖像的曝光時間快速選擇方法進行了研究，本文將針對多曝光圖像的快速融合進行研究。

基于多分辨率分析的像素級融合已經成為圖像融合領域中發(fā)展較快的一種方法［10］。多分辨率分析將圖像分解為能保持局部信息的不同尺度或不同方向的子圖像系列，它們分別代表不同的特征，如邊緣、零交叉、梯度、對比度等，這更能滿足融合的需要?；诙喾直媛史治龅娜诤纤惴ǖ膬?yōu)點主要表現(xiàn)在：能提供對人眼的視覺比較敏感的強對比度信息，以及它在空間和頻域的局部化能力。

基于多分辨率分析的圖像融合方法主要包括兩大類融合方法，即基于塔型變換的融合方法和基于小波變換的融合方法。而其在圖像融合領域中的應用多在多源和多聚焦圖像的融合研究上，對多曝光圖像的快速融合研究較少。多曝光的高動態(tài)范圍視頻系統(tǒng)需要解決以下問題：多曝光圖像的曝光時間快速選擇、快速配準、快速融合，以及偽影去除。

本文著重討論將多曝光應用于HDR視頻快速融合的實時性問題。實時生成的HDR圖像序列可用于直接實時顯示，或作為下一級系統(tǒng)（如視頻壓縮、存儲、或傳輸系統(tǒng)）的輸入。但這些問題不在本文的討論范圍之內。視頻采集合成系統(tǒng)對于計算的實時性有嚴格的要求，如對于常用的25fps系統(tǒng)來說，要求每幀處理時間不大于40ms。

塔形變換分解和合成的計算復雜度比小波變換要小得多。出于對速度的考慮，本文采用金字塔分解的方法對曝光圖像進行融合，并對塔形分解中的重構算法、拉普拉斯系數(shù)融合規(guī)則、高斯系數(shù)融合規(guī)則進行了改進，在不降低融合質量的前提下大大減少了計算量，適應了實時視頻采集。

2 基于拉普拉斯金字塔的快速多曝光融合

2.1 基于拉普拉斯金字塔分解的多曝光融合

拉普拉斯金字塔是Burt提出的一種類似于帶通濾波的空間域技術［12］。算法首先通過一個對稱的低通濾波器w（如5×5）對原始M×N的圖像g0卷積，將得到的圖像進行1／2下采樣，得到一幅（M／2）×（N／2）的圖像g1，這個過程重復L次，得到一系列尺寸逐級減半的圖像g1，g2，...，gL。由于濾波器w的形狀通常類似高斯分布曲線，因此這些分解得到的圖像被稱為高斯金字塔或高斯分解系數(shù)。其計算過程被稱為REDUCE操作，如式（1）所示。

高斯分解系數(shù)可理解為圖像的低頻分量，是圖像能量主要的集中部分。Burt希望通過尺寸較小、但能量很高的低頻分量（最高一級的高斯分量）和一系列能量很低的高頻分量來表示原圖像，從而實現(xiàn)圖像的壓縮編碼。因此他設計了一種計算第i－1級的高頻分量的方法：將第i級高斯分解gi通過擴展插值操作（EXPAND）預測其1／2下采樣前的第i－1級高斯分解gi－1，得到g＿predi－1，并計算其預測誤差，即為第i－1級的拉普拉斯分量lpi－1：

拉普拉斯分量可以理解為圖像的高頻分量，包含了圖像的邊緣紋理細節(jié)信息。

而原始圖像的重建過程只是分解的反向操作。

文獻［13］描述了一種基于拉普拉斯金字塔分解的多曝光融合算法（稱之為Mertens算法）。設有兩個已配準的不同曝光時間的M×N大小的圖像I1，I2，并計算兩幅圖像的權重分布圖W1，W2。將兩幅圖像進行L級拉普拉斯金字塔分解，分別得到L個拉普拉斯分量lpI1（0）～lpI1（L－1），lpI2（0）～lpI2（L－1），和1個最高級的高斯分量gsI1（L），gsI2（L）。再將權重分布圖 W1，W2也進行L級拉普拉斯分解，但只計算其L個高斯分量gsW1（1）～gsW1（L），gsW2（1）～gsW2（L），特殊地，第0級高斯分量就是W1，W2：gsW1（0）＝W1，gsW2（0）＝W2。

其權重函數(shù)采用高斯函數(shù)：

用gsW1（0）～gsW1（L），gsW2（0）～gsW2（L）對lpI1（0）～lpI1（L－1）、gsI1（L）和lpI2（0）～lpI2（L－1）、gsI2（L）的對應等級進行加權融合，如式（5）所示。

其中l(wèi)pF（i）（x，y）代表融合后的第i級拉普拉斯分解在坐標（x，y）處的系數(shù)；gsF（L）（x，y）代表融合后的最高級高斯分解在坐標（x，y）處的系數(shù)；P為待融合的圖像數(shù)量，本例中P＝2。

這樣就得到了融合后的L級拉普拉斯分量lpF（0）～lpF（L－1）和1個最高級高斯分量gsF（L）。再通過重構，得到融合后的圖像。

通過實驗得，一般分解等級L＝3就可以得到較好的融合效果。

2.2 改進的快速塔形分解多曝光融合

上述方法雖然應用了拉普拉斯塔形變換，避免了小波變換帶來的復雜度提高，但其計算量仍然巨大，難以將計算時間減少到40ms內，無法適應實時的視頻采集融合。

本節(jié)通過4個方面的改進，在保證融合質量的前提下，大幅降低運算復雜度，提高計算速度，滿足實時性要求。

2.2.1 金字塔分解卷積核的改進

Mertens應用的拉普拉斯金字塔分解使用的卷積核為5×5的矩陣：

但由于這樣的卷積核在分解塔形結構時將在每個像素上計算25次乘法和24次加法，以及最后除以256的25次移位操作。以512×512的灰度圖進行3級分解為例，將計算的次乘法、加法和移位次數(shù)約為：

若只對兩幅圖像分解，對于4800MIPS的TMS320DM642，不考慮存儲器延時及流水線破壞等因素，僅高斯分解就將耗時2～3ms。

因此考慮降低卷積計算復雜度，提高卷積速度。將卷積核修改為以下3×3的矩陣：

這樣每個像素的卷積將變?yōu)?次乘法、8次加法以及最后除以16的9次移位，其中乘法都可由算術移位來實現(xiàn)，且4個角上的元素不需要計算。這樣共計14次移位和8次加法。對兩幅圖像進行高斯分解將僅耗時0.4ms，速度提升5～7倍。

過小的卷積核將使圖像高頻和低頻分解不徹底，導致高層次的高斯分量仍有大量細節(jié)，使融合不理想，這可以通過提高分解級數(shù)來解決。但實驗表明，這種卷積核在3級分解下并沒有對圖像融合產生太大影響。

2.2.2 簡化高斯系數(shù)擴展插值方法（EXPAND）

高斯系數(shù)的擴展插值預測是影響計算速度的重要因素，因為這個過程在金字塔分解和重構的過程中都會存在。原有的預測方法如式（3）所示。這是個類似卷積的操作，使用與分解過程相同的卷積核。以一維、5抽頭濾波器為例，其插值關系如圖1所示，由第i級高斯系數(shù)插值得到第i－1級預測。

圖1 拉普拉斯金字塔插值操作示意圖Fig.1 Sketch map of interpolating method for Laplacian Pyramid

可以看到每插值一個像素，需要計算5×5＝25次乘法、24次加法和最后除以256的25次移位。以3級分解512×512圖像為例，分解過程中的插值操作需要的乘法、加法和移位操作次數(shù)約為：

而重構需要同樣的計算量。對于4800 MIPS的TMS320DM642來說不考慮存儲器延時及流水線破壞等因素，僅分解和重構的插值操作就將耗時10ms。

圖2 簡單直接插值預測產生的方塊效應Fig.2 Square distortion caused by direct interpolate prediction

因此需要改進插值算法，降低差值算法的復雜度，提高運算速度。最簡單的插值方法就是用當前層次的高斯系數(shù)直接給上一層的對應位置的2×2像素賦值，但這樣做會在紋理邊緣產生嚴重的預測誤差，使融合圖像出現(xiàn)方塊效應，且隨著分解級數(shù)的提高而越來越嚴重。如圖2所示。

圖3 簡化的插值操作示意圖Fig.3 Sketch map of simplified interpolating

我們將插值操作修改為如圖3的形式，第i級高斯分量系數(shù)gsi（0）對應的4個第i－1級高斯分量系數(shù) gsi－1（0），gsi－1（1），gsi－1（2），gsi－1（3），按如下方式計算：

這在一定程度上減小了預測誤差，緩解了方塊的產生，如圖4所示，可以看到，融合圖像中已不存在方塊效應。

圖4 簡單平均值預測融合結果Fig.4 Results of simplified average prediction merge

經過這樣的改進，一個像素的插值預測過程將只有3次加法和3次移位，分解過程中的插值操作共需1.0×106次加法和移位，重構同理。這樣對于4800MIPS的TMS320DM642來說不考慮存儲器延時及流水線破壞等因素，分解和重構中的插值操作將僅耗時約0.8ms，提高了12倍。

實驗表明，簡化的預測方法與原有的方法相比，圖像質量沒有太大的改變。

2.2.3 簡化各級拉普拉斯高頻系數(shù)融合方法

為了將各級拉普拉斯系數(shù)相融合，Mertens將兩幅圖像的權重圖也進行了塔形分解，以得到的各級高斯分量作為待融合圖像的各級拉普拉斯分量和最后一級高斯分量的權重圖，其基本思想是認為像素值在中間位置的像素所產生的紋理信息更為可信，從而賦予其更高的權重。如圖5（a）所示。

但這樣融合2幅多曝光圖像除了需要對2幅原圖像進行塔形分解，還要對其2幅權重圖進行塔形分解，共需4次塔形分解和1次塔形重構。因此我們嘗試簡化各級權重圖的獲取方法：根據待融合原圖像的各級高斯分量直接按照權重函數(shù)查表獲得各自的權重圖，如圖5（b）所示。由于在塔形分解過程中會產生各級高斯分量，因此可以充分利用分解過程中所產生的計算結果，避免對權重圖的塔形分解。

圖5 （a）Mertens計算塔形分解各級權重的方法；（b）權重不進行塔形分解.Fig.5 （a）Mertens’method of weight map computing for each pyramid level；（b）The weight map is not pyramid decomposed.

實驗表明這種權重圖的生成方法與圖5（a）的方法相比，融合結果沒有產生重大差異。

然而實驗中也發(fā)現(xiàn)這兩種方法融合得到的圖像邊緣紋理出現(xiàn)模糊，這是對高頻系數(shù)加權時產生的必然結果。因此考慮在保證邊緣細節(jié)的同時進一步大幅減少計算量。本文對下一種方法進行了實驗：高頻系數(shù)融合時只選擇絕對值最大的系數(shù)，以充分保留圖像的細節(jié)部分，并減少計算量。

實驗結果表明效果明顯，得到的圖像邊緣紋理清晰，相比Mertens算法，圖像質量有了較大提高。但這種方法也有缺點，就是對噪聲比較敏感，要求輸入圖像源噪聲比較低。

2.2.4 簡化高斯低頻系數(shù)融合方法

Mertens在對最低分辨率的低頻系數(shù)融合時使用了分解得到的權重圖。與2.3.3所敘述的一樣，這種方法存在多次塔形分解，耗費計算量的問題。應用2.3.3中所提到的第一種方法，計算量得到了明顯降低，但實驗過程中也發(fā)現(xiàn)低頻的融合對權重函數(shù)非常依賴，且容易引起亮度分布的不均，出現(xiàn)斑塊。

我們直接將低頻系數(shù)取平均值，使亮度分布均勻，避免出現(xiàn)斑塊，同時也大大降低了計算量。實驗表明亮度分布不均的問題得到了解決，但高輻照度區(qū)域的對比度降低了，這有待進一步改進。

2.2.5 改進后的算法計算量分析

通過以上一系列的改進，多曝光的融合速度得到了很大的提高，而融合質量并沒有因為算法的簡化而發(fā)生降低，反而在圖像細節(jié)融合質量上有些提高。但高輻照度區(qū)域的對比度有些下降。

以512×512大小的圖像在TMS320DM642上實現(xiàn)快速融合為例，我們對本文快速塔形分解算法和Mertens算法的計算量做了粗略統(tǒng)計，如表1所示。這個統(tǒng)計沒有考慮流水線的阻塞、存儲器延時及輔助指令等因素，只統(tǒng)計了主要的運算指令，但這足以說明計算量的大致關系。

從表1中可以看到，Mertens算法中較為耗時的操作為Reduce、Expand和加權平均操作。

本文重點改進了以上3種操作，使其計算量大為降低，總計算量削減了90%，計算速度提高了10倍。

最后等效計算時間按照未進行并行優(yōu)化的代碼比較，Mertens算法主要運算指令時間為192 ms，即使按照100%并行（每個指令周期執(zhí)行8條指令，通常無法實現(xiàn)）折算，實際硬件計算時間也要達到24ms，若考慮流水線的阻塞、存儲器延時及輔助指令等因素，無法實現(xiàn)40ms的實時視頻合成計算要求。而這些估算還僅僅是對于512×512的灰度視頻圖像序列，對于更大尺寸視頻序列，Mertens算法更不能滿足其融合的實時性要求。

表1 Mertens算法與本文快速塔形分解算法計算量比較（單位512×512cycle）Tab.1 Calculation amount comparison between Mertens and proposed LDR image fusion method （unit：512×512cycle）

3 實驗結果與分析

3.1 系統(tǒng)硬件結構

本文算法被應用于某運載火箭視頻采集系統(tǒng)，其系統(tǒng)框圖和實物圖6所示。

視頻捕捉部分為BASLAR CCD數(shù)字相機，輸出LVDS信號，通過串并轉換器轉換為并行數(shù)據輸入給DM642的視頻接口VP0。同時DM642通過RS232實現(xiàn)CCD數(shù)字相機的快速曝光時間調節(jié)，獲得多曝光視頻序列。得到多曝光圖像數(shù)據后由DM642執(zhí)行本文快速融合算法，并通過視頻口VP1將融合后的數(shù)據發(fā)送給視頻壓縮及傳輸系統(tǒng)，將高動態(tài)范圍視頻傳輸至地面測控系統(tǒng)。

圖6 某運載火箭高動態(tài)范圍視頻采集融合系統(tǒng)框圖及實物Fig.6 System diagram and picture of a launch vehicle HDR video observation system

為方便地面實驗及算法驗證，融合后的數(shù)據還通過以太網將數(shù)據發(fā)送給PC機進行實時顯示（虛線部分）。

3.2 多曝光融合的質量評價實驗

圖7 各種融合方法對Garden的融合結果Fig.7 Fusion results of Garden using different methods

為了驗證本文融合方法的有效及快速性，分別對本文中提到的兩種方法進行對比，所采用的視頻序列使用曝光時間Δt＝0.15和Δt＝2.20的Garden（493×874）、曝光時間 Δt＝1／2090s和Δt＝1／159s的真實拍攝場景 EsteeLauder（550×680）、Δt＝1／60s和Δt＝1／6000s的真實拍攝場景DongPei（550×680）。圖7是Garden的融合結果。

可以看到，輻照圖重建法不能達到較好的融合效果。而加權融合丟失了很多細節(jié)。

圖8是EsteeLauder的融合結果。

圖8 各種融合方法對EsteeLauder的融合結果Fig.8 Fusion results of EsteeLauder using different methods

圖9 各種融合方法對Dongpei的融合結果Fig.9 Fusion results of Dongpei using different methods

圖9是對Dongpei的融合結果。

在主觀質量上，以上實驗結果表明本文快速塔形分解算法融合的圖像比其他算法都清晰銳利。

為了對融合結果有一個客觀的評價，引入了圖像的信息熵。圖像的信息熵表征了像素值在取值范圍內的離散程度。過曝和欠曝的圖像都會引起圖像直方圖的不均勻分布，從而導致信息熵的低下。而紋理、邊緣細節(jié)豐富的圖像趨于信息熵的最大化。但圖像的大量噪聲也會引起信息熵的提高，不能一概用熵來評價。但可用其作為一個參考，輔以人的主觀感受，對圖像質量有一個合理的評價。

表2是對各拍攝場景不同融合結果的信息熵客觀評價。其中列向灰色底紋的單元格表示對這一場景下的熵較高的兩種融合方法。

可以看到快速塔形分解算法在多個測試實驗中都能夠保持較高的信息熵，且都要高出Mertens算法，而其他方法測得的熵并不是很穩(wěn)定，只對某一種場景可以得到較高的熵，因此本文的方法穩(wěn)定性更強，且無論在主觀感受和客觀評價上融合效果都比Mertens算法略好。

表2 不同場景使用不同的融合方法的信息熵評價Tab.2 Entropy appraising of different fusion method for different scene

3.3 TMS320DM642實現(xiàn)的計算時間比較

我們將Mertens算法和本文快速塔形分解算法在TMS320DM642上進行了計算，利用CCS的Profile功能對程序執(zhí)行情況作了統(tǒng)計，結果如表3所示。

表3 快速塔形分解算法和Mertens算法在TMS320DM642＠600MHz的計算時間Tab.3 Computing time on TMS320DM642＠600MHz of Mertens’and proposed method

其中，CPU 總阻塞周期（Stall Cycles）包含L1DCache阻塞周期（L1DStall Cycles）。

可以看到總計算時間縮短到了8.9ms，提高了10.6倍，這是由于改進了塔形分解和重構的REDUCE、EXPAND算法、去除了權重圖的塔形分解等造成的。

其中Mertens法由于L1Cache缺失引起的流水阻塞周期（L1DStall Cycles）比本文快速塔形分解算法多很多，這是由于5×5的卷積和使每計算一個高斯系數(shù)需要讀取更多的像素，產生更多的讀存儲器指令。

同時也可以看到核心計算時間（除去阻塞周期的計算周期，Core Cycles）也約為快速塔形分解算法方法的8.6倍，與預計結果稍有出入，這是由于在計算中不單有對像素的計算，還包含了存儲單元的存取、索引等控制代碼，以及單指令多數(shù)據（SIMD）指令給統(tǒng)計結果帶來的不確定影響。

4 結 論

在使用多曝光法進行高動態(tài)范圍視頻的合成時，計算的實時性是急需解決的問題。Mertens提出了一種基于拉普拉斯金字塔分解的多曝光融合算法，但其無法滿足視頻采集融合系統(tǒng)的實時性要求。

本文對塔形分解中的重構算法、拉普拉斯系數(shù)融合規(guī)則、高斯系數(shù)融合規(guī)則進行了改進，在不明顯降低融合質量的前提下盡量減少計算量，以適應實時多曝光視頻融合需求。

通過硬件系統(tǒng)實驗，驗證了該改進算法的穩(wěn)定性及融合質量，同時減少了計算量。對于512×512的灰度視頻序列，運算時間在600MHz的C6000系列DSP上由94ms減少到了8.9ms，計算速度提高了10倍，滿足了一般視頻采集融合系統(tǒng)25fps的實時性要求。

［1］關澈，王延杰.CCD 相機實時自動調光系統(tǒng)［J］.光學精密工程，2008，16（2）：358－366.Guan C，Wang Y.Real－time auto light control system of CCD camera［J］.Opt.Precision Eng.，2008，16（2）：358－366.（in Chinese）

［2］楊海濤，常義林，王靜，等.一種基于亮度直方圖的自動曝光控制方法［J］.光學學報，2007，27（5）：841－847.Zhou H T，Chang Y L，Wang J，et al.A new automatic exposure algorithm for video cameras using luminance histogram ［J］.Acta Optica Sinica，2001，8（6）：50－53.（in Chinese）

［3］蘇宏武，楊小君，李哲，等.基于平均和峰值灰度加權的自動調光系統(tǒng)［J］.光子學報，2006，35（1）：158－160.Su H W，Yang X J，Li Z，et al.Automatic adjust light systems based on combined average gray with peak gray［J］.Acta Optica Sinica，2006，35（1）：158－160.（in Chinese）

［4］Mann S，Picard R W.On being‘undigital’with digital cameras：extending dynamic range by combining differently exposed pictures［C］.IS＆T’s 48th Annual Conference，1995：422－428.

［5］Nayar S K，Mitsunaga T.High dynamic range imaging：spatially varying pixel exposures［C］.Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition，2000.

［6］Robertson M A，Borman S，Stevenson R L.Estimation－theoretic approach to dynamic range enhancement using multiple exposures［J］.Journal of Electronic Imaging，2003，12（2）：219－228.

［7］楊作廷，阮萍，翟波.基于圖像熵的高動態(tài)范圍場景的自動曝光算法［J］.光子學報，2013，42（6）：742－746.Yang Z T，Ruan P，Zhai B.Auto－exposure algorithm for scenes with high dynamic range based on image entropy［J］.Acta Optica Sinica，2013，42（6）：742－746.（in Chinese）

［8］唐艷秋，張星祥，李新娥，等.基于人眼視覺灰度識別特性的圖像動態(tài)范圍小波變換處理方法［J］.液晶與顯示，2012，27（3）：385－390.Tang Y Q，Zhang X X，Li X E，et al.Image processing method of dynamic range with wavelet transform based on human visual gray recognition characteristics［J］.Chinese Journal of Liquid Crystals and Displays，2012，27（3）：385－390.（in Chinese）

［9］樸永杰，徐偉，張貴祥.高動態(tài)范圍視頻的多曝光時間自適應快速選擇［J］.光學精密工程，2012，20（12）：216－223.Piao Y J，Xu W，Zhang G X.Fast and adaptive decision of exposure time in multi－exposure HDR video［J］.Opt.Precision Eng.，2012，20（12）：216－223.（in Chinese）

［10］Piella G.A region－based multiresolution image fusion algorithm ［C］.Proeeedings of the 5thInternational Conference on Information Fusion，2002：1557－1564.

［11］Barakat N，Darcie T E.Minimal capture sets for multi－exposure enhanced－dynamic－range imaging ［C］.Proceedings of International Symposium on Signal Processing and Information Technology，2006：524－529.

［12］Burt P J，Anelson E H.The Laplacian pyramid as a compact image code［J］.IEEE Transactions on Communications，1983，31（4）：532－540.

［13］Mertens T，Kautz J，Reeth F V.Exposure fusion ［C］.Proceedings of 15th Pacific Conference on Computer Graphics and Applications，2007：382－390.

［14］Schlik C.Quantization techniques for visualization of high dynamic range pictures［C］.5th Eurographics Workshop on Rendering，1994：7－18.

［15］Vuong Q K，Yun S H，Kim S.A new auto exposure and auto white－balance algorithm to detect high dynamic range conditions using CMOS technology［C］.Proceedings of the World Congress on Engineering and Computer Science，2008：210－205.

［16］Blinn J.Dirty pixels［J］.IEEE Computer Graphics Application，1989（9）：100－105.