劉凱，寇正，羅立民

1.東南大學(xué)計算機科學(xué)與工程學(xué)院，南京 210096

2.解放軍理工大學(xué)氣象海洋學(xué)院，南京 211101

PDE融合方法的兩種離散方案的比較研究

劉凱1，2，寇正2，羅立民1

1.東南大學(xué)計算機科學(xué)與工程學(xué)院，南京 210096

2.解放軍理工大學(xué)氣象海洋學(xué)院，南京 211101

1 引言

近年來，圖像融合技術(shù)特別是像素級圖像融合技術(shù)得到了迅速的發(fā)展，從經(jīng)典的基于金字塔變換的方法[1-4]到基于多尺度幾何分析的方法[5-8]，可以說圖像融合算法層出不窮。Pop.S等[9-11]從2007年開始相繼發(fā)表論文獨創(chuàng)性地提出了基于偏微分方程的圖像融合方法，他們認(rèn)為圖像融合過程就是對每一幅源圖像應(yīng)用一個基于偏微分方程的演化過程，在演化過程的每一步，保留當(dāng)前源圖像中的相關(guān)信息的同時加入其他源圖像中的信息，為此他們提出了描述每一幅源圖像演化過程的一般連續(xù)方程。在與主成分分析（PCA）方法和多分辨率分析（MRA）方法等經(jīng)典融合方法進(jìn)行對比實驗的結(jié)果顯示，基于偏微分方程的融合方法生成的融合圖像的質(zhì)量或比經(jīng)典方法高或與經(jīng)典方法接近。他們的工作使得基于偏微分方程的圖像融合方法成為一類新興的融合方法，并開辟了偏微分方程應(yīng)用的新領(lǐng)域，值得廣大科研工作者進(jìn)行進(jìn)一步深入的探討和研究。在連續(xù)方程被應(yīng)用之前必須對其進(jìn)行離散化，Pop.S等采取的離散化方案是逆擴(kuò)散項采用先計算梯度再求散度的方案，而本文就在原始連續(xù)方程的基礎(chǔ)上，通過變換和對方程的正則化項進(jìn)行相應(yīng)的調(diào)整，推導(dǎo)出另一種形式的方程，將原先的離散化方案變化為計算拉普拉斯的離散化方案，并且將兩種方案進(jìn)行比較。

2 基于偏微分方程的圖像融合算法

Pop.S等提出了一個一般連續(xù)的源圖像的演化方程。

2.1 一維連續(xù)方程

方程（1）右邊第一項稱為融合項，也被稱為逆擴(kuò)散項，其目的是將其他源圖像中的相關(guān)信息加入到當(dāng)前源圖像中。式中，div（）是散度算子，i代表當(dāng)前源圖像，max代表具有最大梯度絕對值的源圖像，βi為正的權(quán)系數(shù)。

對當(dāng)前像素而言，如果最大梯度絕對值出現(xiàn)在當(dāng)前源圖像中，則當(dāng)前像素不變(βi=0)；否則，如果最大梯度絕對值在其他源圖像中被檢測到，則梯度以逆擴(kuò)散過程的方式加入到當(dāng)前圖像中。

方程（1）右邊第二項為正則化項，其目的是避免模型產(chǎn)生振蕩，因為如果方程中只有逆擴(kuò)散項可能會導(dǎo)致振蕩和不穩(wěn)定，需要對方程施加約束條件，約束條件就以正則化項表示。式中，γ是正的正則化權(quán)系數(shù)，其設(shè)定正則化項的重要性，gR是一個分段線性函數(shù)。正則化項的形式完全取決于融合項，若令gR=1，其形式完全與融合項等同，它的作用相當(dāng)于一個梯度極值約束（文獻(xiàn)[9-11]），使積分過程中方程右端的量值不超過完全由初始梯度極值計算的融合項。

同時方程（1）也表示了一幅圖像的演化過程，此方法的機制是每一幅輸入的源圖像都進(jìn)行這樣的演化過程，在經(jīng)過一個時間步長后所有圖像同步更新數(shù)據(jù)。區(qū)別與傳統(tǒng)融合方法的是此算法使得每個源信號都會產(chǎn)生一個輸出，當(dāng)所有輸出幾乎達(dá)到相等時，可以認(rèn)為融合完成。因此融合的過程實際上也是輸出信號之間均方根誤差（RMSE）隨時間不斷減小的收斂過程。

2.2 數(shù)值化方案

方程（1）的離散化方案為：

ma代表當(dāng)前像素的最大前向和最大后向差分絕對值所在的源圖像。

方程（3）顯示了div(?Umax)項的離散化過程，首先，對當(dāng)前像素檢測最大前差絕對值(ma)和最大后差絕對值(ma)，如果ma或ma發(fā)生在當(dāng)前源圖像中，則=0或=0；否則，∈[0,1]或∈[0,1]?？紤]到在所研究

的鄰域內(nèi)，ma和ma可能出現(xiàn)在不同的源圖像中，那么，βi就可能有兩個值。然后，前項減后項得到散度。

3 新的融合方法

3.1 連續(xù)方程

根據(jù)數(shù)學(xué)定義，方程（1）的右邊第一項就是：

式中，?2()是拉普拉斯算子，max代表具有最大拉普拉斯絕對值的源圖像。

方程的正則化項也根據(jù)第一項進(jìn)行相應(yīng)調(diào)整，這樣方程（1）就變換為：

方程（9）的第一項和第二項仍然為融合項和正則化項，但與方程（1）不同的是，方程（9）反映了融合圖像中的每一像素的信息是由具有最大拉普拉斯絕對值的源圖像提供。參照方程（1）中第二項的形式，也給方程（9）添加一個拉普拉斯極值約束，即方程（9）右端第二項。

3.2 數(shù)值化方案

[12-13]中的離散化方案，拉普拉斯算子?2Umax的離散為：

方程（12）說明了在新方法中對div(?Umax)的離散化過程，就是對當(dāng)前像素直接檢測拉普拉斯的最大絕對值()，如果出現(xiàn)在當(dāng)前源圖像中，βi=0；否則，βi∈[0,1]。

兩種融合方法的區(qū)別就在于對逆擴(kuò)散項的離散化方案，原方法中是通過先計算梯度再求解散度的方案，而在新方法中是直接計算拉普拉斯。

數(shù)值積分過程均采用了固定邊界條件。

4 兩種方案的比較與分析

4.1 一維情形下的比較

為了清楚地反映兩種方案的不同作用效果，將兩種方案分別應(yīng)用于兩個一維信號的融合，為什么選擇在一維情形下進(jìn)行比較，因為一維信號簡單、只涉及一個方向，可以清楚地看到兩種方案在運行過程中的差異。這里假定信號A為零，信號B分別為方波、梯形波和三角波，使用兩種方案將信號A和信號B進(jìn)行融合。圖1和圖2分別是應(yīng)用原方案和新方案將零與方波融合時信號的演化過程。圖3是兩種方案在信號融合過程中兩信號均方根誤差的變化曲線。圖4和圖5分別是應(yīng)用原方案和新方案將零與梯形波融合時信號的演化過程。圖6是兩種方案在信號融合過程中兩信號均方根誤差的變化曲線。圖7和圖8分別是應(yīng)用原方案和新方案將零與三角波融合時信號的演化過程。圖9是兩種方案在信號融合過程中兩信號均方根誤差的變化曲線。上述所有一維信號實驗的參數(shù)設(shè)置是相同的，即gF=1,β=1,γ=1,dt=0.1。

圖1 應(yīng)用原方案將信號A與方波融合的信號演化過程

2 應(yīng)用新方案將信號A與方波融合的信號演化過程

圖3 兩方案對零和方波信號融合時的均方根誤差變化曲線

圖4 應(yīng)用原方案將信號A與梯形波融合的信號演化過程

圖6 兩方案對零和梯形波信號融合時的均方根誤差變化曲線

由圖1和圖2所示，兩種方案對零與方波的融合演化過程中在各時刻信號的變化和波形都是十分相似的，當(dāng)t= 30時，基本都完成了兩信號的融合。再由圖3所示，兩方案的RMSE變化曲線幾乎重合，收斂速度幾乎相等，可以說兩種離散化方案對零與方波的融合幾乎沒有差別。但當(dāng)信號B為梯形波時，由圖4和圖5所示，應(yīng)用新方案比原方案融合的速度快。對新方案，當(dāng)t=30時，信號A和信號B都演化成梯形波的波形；而與此同時，對原方案，信號A和信號B都還沒有演化到梯形波的波形。這一點也由圖6得到證實，如圖6所示，新方案的RMSE比原方案下降得快，在t=30時已經(jīng)下降到0，新方案加速了收斂過程。當(dāng)信號B為三角波時，情形與梯形波時類似，由圖7～圖9所示，應(yīng)用新方案也是比原方案融合的速度快，新方案明顯加速了收斂過程。

圖5 應(yīng)用新方案將信號A與梯形波融合的信號演化過程

從上面的三組信號的實驗結(jié)果可以發(fā)現(xiàn)，兩種離散化方案對突變信號（例如方波）的融合效果基本相同，差別很小。但對緩變信號（例如梯形波、三角波）時，新方案比原方案收斂速度要快。

4.2 二維情形下的比較

將兩種方案擴(kuò)展到二維，分別從x和y兩個方向進(jìn)行計算，比較兩種方案對二維圖像的融合效果。選擇了兩組經(jīng)常用于融合實驗的源圖像，一組是聚焦不良圖像，另一組是醫(yī)學(xué)圖像。參考文獻(xiàn)[11]將新方案與原方案從融合圖像質(zhì)量指標(biāo)和算法的收斂性兩個方面進(jìn)行比較，為了定量比較融合圖像的質(zhì)量，選用Piella等[14]提出的加權(quán)融合圖像質(zhì)量評價指標(biāo)Qw，Qw值越大，融合圖像的質(zhì)量越高。算法的收斂性比較上，仍然采用輸出圖像的均方根誤差（RMSE）的演變曲線來分析、比較，RMSE下降越快，說明算法的收斂越快；RMSE值越小，說明輸出的圖像越相似。

圖7 應(yīng)用原方案將信號A與三角波融合的信號演化過程

圖8 應(yīng)用新方案將信號A與三角波融合的信號演化過程

圖9 兩方案對零和三角波信號融合時的均方根誤差變化曲線

第一組實驗的參數(shù)設(shè)置：β=0.1，γ=2.4,dt=0.1。第二組實驗的參數(shù)設(shè)置：β=0.5，γ=2,dt=0.05。

圖10（a）、圖10（b）和圖11（a）、圖11（b）分別為第一組和第二組融合實驗源圖像。圖10（c）、圖10（d）和圖11（c）、圖11（d）分別是采用原方案和新方案經(jīng)過1 000次迭代后得到的融合圖像。圖10（e）、圖10（f）和圖11（e）、圖11（f）分別是兩組源圖像的兩種方法的Qw隨迭代次數(shù)N和RMSE隨迭代次數(shù)N的變化曲線，其中實線代表原方案，虛線代表新方案。

圖10 第一例融合實驗結(jié)果

圖11 第二例融合實驗結(jié)果

就第一例實驗而言，由圖10（c）和圖10（d）所示，相對于源圖像，無論是原方案還是新方案經(jīng)過1 000次迭代后輸出的融合圖像都取得了良好的視覺效果，源圖像中模糊的區(qū)域經(jīng)過融合后都消除了模糊，變得清晰了。從質(zhì)量指標(biāo)Qw看，由圖10（e）中的實線和虛線所示，剛開始隨著迭代次數(shù)N的增加，兩種方案的Qw都迅速上升，提高很快，而新方案的虛線更陡直，即新方案的Qw上升速度更快一些，Qw值也略高，說明采用新方案融合圖像的質(zhì)量比采用原方案的圖像質(zhì)量略好，同時經(jīng)過較少的迭代次數(shù)，就能達(dá)到和采用原方案同樣的質(zhì)量指標(biāo)。當(dāng)?shù)螖?shù)N超過300后，隨著迭代次數(shù)N的增加，Qw都上升緩慢，并且在此之后近似保持不變，但采用新方案的Qw值也總是比采用原方案的略高。而且從兩曲線的走勢看，隨著迭代次數(shù)的繼續(xù)增加超過1 000次后，可以推斷采用新方案的融合圖像的質(zhì)量始終比采用原方法的圖像質(zhì)量略好。再從算法的收斂性上看，由圖10（f）中的實線和虛線所示，剛開始隨著迭代次數(shù)N的增加，兩種方案的RMSE都迅速下降，但虛線比實線下降得更快，RMSE值更小，說明采用新方案輸出圖像收斂更快，輸出圖像更相似；同樣當(dāng)?shù)螖?shù)N超過300后，兩種方案的RMSE都下降非常緩慢，并基本保持不變，但采用新方案的RMSE值總是比采用原方案的RMSE值小，說明如果經(jīng)過相同次數(shù)的迭代，采用新方案的輸出圖像比采用原方案的輸出圖像更相似。

對第二例實驗而言，由圖11（c）和圖11（d）所示，同樣相對于源圖像，兩種方案下經(jīng)過1 000次迭代后輸出的融合圖像也都取得了良好的視覺效果。由圖11（e）中實線和虛線所示，經(jīng)過1 000次迭代后，采用兩種方案的融合圖像的Qw都很接近于最大值1。但是在這之前，虛線比實線更陡直一些，采用新方案的Qw上升得略微快一些，說明采用新方案經(jīng)過較少的迭代次數(shù)N，就能達(dá)到較高的質(zhì)量指標(biāo)。再從算法的收斂性上看，經(jīng)過1 000次迭代后，采用新方案的RMSE隨著迭代次數(shù)N的增加，下降得更快一些，輸出圖像也更相似，說明其收斂性更好。

從兩例實驗的結(jié)果看，總的來說，在相同的參數(shù)設(shè)置條件下，不論從質(zhì)量指標(biāo)還是收斂性上，新方案都比原方案有所提高。

4.3 兩種方案得到不同結(jié)果的原因分析

為什么對某些信號新方案比原方案結(jié)果好?就融合項和正則化項分別進(jìn)行分析。首先看融合項，在連續(xù)情況下，-?2U=-div(?U)，而在原方案中相應(yīng)融合項離散化形式為：，新方案離散化形式為：發(fā)生在不同源圖像中，原方案前后兩項則包含了不同源圖像的信息，也就不嚴(yán)格等于拉普拉斯了。然而在新方案中，直接檢測當(dāng)前像素鄰域內(nèi)拉普拉斯的絕對值()，該項必然只包含一個源圖像的信息，是嚴(yán)格意義的拉普拉斯項。再看正則化項，正則化項前面已說明，它的作用是梯度（拉普拉斯）極值約束。如果在原方案中保持融合項不變，改變正則化項，采用Laplace正則化，以信號A與梯形波融合為例，如圖12所示，兩信號也可以很好地融合，但其RMSE變化曲線與原方案非常接近，沒有明顯改進(jìn)。而在新方案中融合項采用Laplace，正則化項采用Laplace正則化，如圖12所示，新方案的RMSE明顯比原方案下降，這說明融合項的新的離散化方案是導(dǎo)致新方案比原方案優(yōu)越的主要原因，而兩方案中的正則化項雖然形式不同，但不是造成兩種方案差異的根本原因。

圖12 均方根誤差變化曲線

5 結(jié)束語

在本文中，通過對原始連續(xù)方程的變換，提出了基于偏微分方程的圖像融合新的離散化方案，新方法的融合策略是融合圖像的每一像素的信息，是由處于該位置的具有最大拉普拉斯絕對值的源圖像的像素提供的。在一維信號和二維信號融合實驗的結(jié)果顯示，在參數(shù)設(shè)置相同條件下，新方案一般都比原方案的收斂要快，輸出圖像的質(zhì)量指標(biāo)也略高。

參考文獻(xiàn)：

[1]Burt P，Adelson E.Laplacian pyramid as a compact image code[J].IEEE Transactions on Communications，1983，31（4）：532-549.

[2]Burt P.The pyramid as a structure for efficient computation[J]. Multi-resolution Image Processing and Analysis，1984，12.

[3]Toet A，Ruyven L J，Valeton J M.Merging thermal and visual image by a contrast pyramid[J].Optical Engineering，1989，28（7）：789-792.

[4]Toet A.Image fusion by a ratio of low-pass pyramid[J].Pattern Recognition Letters，1989，9（4）：245-253.

[5]Candès E J，Demanet L，Donoho D L，et al.Fast discrete curvelet transforms[J].Multiscale Modeling&Simulation，2005，5（3）：861-899.

[6]Do M N，Vetterli M.The contourlet transform：an efficient directional multi-resolution image representation[J].IEEE Transactions on Image Processing，2005，14（12）：2091-2106.

[7]da Cunha A L，Zhou Jianping，Do M N.The nonsubsampled contourlet transform：theory，design，and applications[J].IEEE Transactions on Image Processing，2006，15（10）：3089-3101.

[8]賈建，焦李成，孫強.基于非下采樣Contourlet變換的多傳感器圖像融合[J].電子學(xué)報，2007，35（10）：1934-1938.

[9]Pop S，Terebes R，Borda M，et al.A PDE-based approach for image fusion[C]//Lecture Notes in Computer Science，2007：121-131.

[10]Pop S，Lavialle O，Donias M，et al.A PDE-based approach to 3D seismic data fusion[J].IEEE Transactions on Geoscience and Remote Sensing，2008，46（5）：1385-1393.

[11]Pop S，Terebes R，Borda M，et al.A new discrete PDE-based fusionmodel[C]//EURCON2009，EURCON’09，2009：2075-2082.

[12]陳利霞，宋國鄉(xiāng)，丁宣浩，等.改進(jìn)的總變分去噪算法[J].光子學(xué)報，2009，38（4）：1001-1004.

[13]周千，趙芳玲，趙鳳群，等.You-Kaveh圖像去噪模型擴(kuò)散系數(shù)的改進(jìn)[J].計算機工程與科學(xué)，2010，32（2）：85-87.

[14]Piella G.New quality metric for image fusion[C]//The 7th International Conference on Information Fusion，Stockholm，Sweden，2004：542-546.

LIU Kai1，2,KOU Zheng2,LUO Limin1

1.School of Computer Science&Engineering,Southeast University,Nanjing 210096,China
2.Institute of Meteorology,PLA University of Science and Technology,Nanjing 211101,China

On the basis of the continuous equation of image fusion based on Partial Differentiation Equations（PDE）proposed by Pop.S et al,a new numerical scheme of the equation is presented.The former discrete scheme of divergence by the calculation of gradient is changed to the scheme of the calculation of Laplace,and the regularization term of equation is also adjusted accordingly.The comparison experiments are designed for two schemes on some 1D and 2D signals.The results show that the new scheme can perform the fusion effectively for all the tested signals;and the new scheme can accelerate the convergence process and make the quality of fused image improved by means of a quality factor for some of the tested 1D and 2D signals.However for the other tested signals,the convergence and the quality of output signals obtained by both the new scheme and the former scheme are also very much closed.

image fusion;partial differential equation;discrete scheme;Laplace;gradient

基于Pop.S等提出的基于偏微分方程（PDE）的連續(xù)圖像融合方程，提出了一種新的離散化方案，將原來的先計算梯度再求散度的離散化方案變化為直接計算拉普拉斯的離散化方案，對方程的正則化項進(jìn)行相應(yīng)的調(diào)整。將這兩種離散化方案分別對一些一維信號和二維信號進(jìn)行融合對比實驗，實驗結(jié)果顯示：對所有的測試信號，新方案都能進(jìn)行有效的融合；對某些測試的一維和二維信號，新方案能夠加速信號的收斂過程并且使得到的融合圖像在質(zhì)量指標(biāo)上比原方案有提高；而對某些測試信號，新方案與原方案在輸出信號的收斂速度上和輸出圖像的質(zhì)量指標(biāo)上也是十分接近的。

圖像融合；偏微分方程；離散化方案；拉普拉斯；梯度

A

TP391.41

10.3778/j.issn.1002-8331.1202-0588

LIU Kai,KOU Zheng,LUO Limin.Comparison study of two numerical schemes for PDE-based fusion method.Computer Engineering and Applications,2013,49（24）：165-171.

劉凱（1970—），女，講師，主要研究方向為氣象遙感衛(wèi)星信息處理、模式識別；寇正（1971—），男，博士，副教授，主要研究方向為大氣動力學(xué)數(shù)值模擬；羅立民，男，教授，博士生導(dǎo)師，主要研究方向為圖像圖形信息處理。E-mail：liukaiqy@163.com

2012-02-29

2012-07-13

1002-8331（2013）24-0165-07