圖像處理中拉普拉斯矩陣的稀疏化處理*

田俊杰，趙祖燁，張軍飛

（1.華中科技大學(xué)材料成型與模具技術(shù)國家重點實驗室，湖北武漢 430074；2.廣州中望龍騰軟件股份有限公司，廣東廣州 510000）

0 引言

圖像處理和計算機(jī)圖形學(xué)中有大量問題可以使用離散泊松方程來求解。圖像處理領(lǐng)域的研究包括Levin等[1]提出的圖像著色，Lischinski等[2]使用的色調(diào)調(diào)整，以及邊緣保留平滑[3]和Xu等[4]提出的相對總變化量紋理剔除算法。泊松方程方法雖然在質(zhì)量和數(shù)學(xué)簡潔性方面表現(xiàn)出色，但其計算成本相當(dāng)較大，需要求解非常大且不易求解的線性系統(tǒng)。

Matrix iterative analysis講述矩陣迭代器運算的基本原理[5]，Iterative methods for sparse linear systems講述求解稀疏線性系統(tǒng)的各種迭代器方法。對于常見圖像處理領(lǐng)域算法中不均勻拉普拉斯矩陣求逆的加速有指導(dǎo)意義[6]。

許多圖像處理算法實現(xiàn)過程中用到稀疏的不均勻拉普拉斯矩陣，當(dāng)圖案像素行列數(shù)分別為M和N時，保邊濾波和圖像著色等算法中的拉普拉斯矩陣的大小為MN×MN。算法求解過程需對拉普拉斯矩陣進(jìn)行求逆，直接求逆的時間復(fù)雜度為O((MN)3)，在時間和內(nèi)存上的消耗都是無法接受的。對于大型稀疏矩陣線性系統(tǒng)可以利用迭代法來加快求解，經(jīng)典的迭代法有Jacobi迭代器方法，Gauss-Seidel迭代器方法，SOR迭代器方法等。本文作者主要應(yīng)用不完全Cholesky分解與預(yù)處理共軛梯度法[6]，并將使用分層稀疏化對迭代器進(jìn)行處理以減少時間和內(nèi)存消耗。

D Krishnan等[7]為在計算機(jī)圖形領(lǐng)域的算法中經(jīng)常出現(xiàn)的離散泊松方程提出一個新的多級稀疏化方案。該方案根據(jù)鄰域內(nèi)粗細(xì)變量的拓?fù)潢P(guān)系選擇不同的稀疏方法，時間復(fù)雜度高，且在分層稀疏化處理過程中無法抑制條件數(shù)的增長，本文作者采用一種更簡單有效的稀疏化方法。

1 算法應(yīng)用場景

Farbman等[3]提出加權(quán)最小二乘濾波算法(WLS)：

式中：g是輸入圖像，找出使式（1）有最小值的u即是想要的輸出圖像；p是像素點的索引；是在像素點p處的梯度；ax,p(u)表示權(quán)重大小。

原理是輸出圖像和輸入圖像的每個像素值盡量接近，且輸出圖像中的非邊界區(qū)梯度盡量小。

式（1）可以轉(zhuǎn)化為矩陣形式：

其中L為大型稀疏的拉普拉斯矩陣，稀疏矩陣求逆一般采用迭代法，但迭代法直接求逆會消耗大量時間。

文獻(xiàn)[4]在WLS的基礎(chǔ)上提出一種相對總變化量算子在保留圖像邊界的同時剔除紋理。紋理和主結(jié)構(gòu)在相對總變化量算子上展現(xiàn)出完全不同的屬性，由此可以在加權(quán)最小二乘法的過程中應(yīng)用不同的權(quán)值對紋理和主結(jié)構(gòu)進(jìn)行不同程度的懲罰。文獻(xiàn)[4]主要的公式是：

其中L也是一個大型稀疏的拉普拉斯矩陣。

2 算法分析

采用一種分層稀疏化的方案減少矩陣求逆的時間消耗。一個矩陣L關(guān)于向量x的能量函數(shù)由Rayleigh商定義:

由于L的半正定型，能量值總是非負(fù)數(shù)。矩陣L的特征向量對應(yīng)的能量值是特征向量的對應(yīng)的特征值，由Lx=λx得(xTLx)/(xTx)=(xTλx)/(xTx)=λ對稱正定矩陣的條件數(shù)定義為：

迭代法通過一次次迭代來逐步逼近離散泊松方程的真實解，而所需的迭代次數(shù)與矩陣的條件數(shù)相關(guān)。Jacobi和Gauss-Seidel迭代法需O(κ)次，共軛梯度法需要O(κ)次[6]。所以加快迭代求解的方式之一就是減小矩陣的條件數(shù)。

關(guān)于減少矩陣條件數(shù)的途徑，D Krishnan等[7]提到一種分層稀疏化方法，每層中都將節(jié)點劃分為粗節(jié)點C和細(xì)節(jié)點F，通過迭代過程分層稀疏化矩陣。矩陣L可以劃分為細(xì)節(jié)點C之間連接LCC，粗節(jié)點F之間連接LFF、細(xì)節(jié)點C與粗節(jié)點F之間連接LFC：

在L兩邊分別乘上得到原問題的一個更小規(guī)模的子問題。

根據(jù)粗細(xì)節(jié)點間的拓?fù)潢P(guān)系選擇不同補(bǔ)償方式稀疏化矩陣，采用一種更簡單有效的矩陣稀疏化方式。采取去除局部三角形中最弱連接，補(bǔ)償給其他兩條邊的方式來稀疏化矩陣以減少條件數(shù)，即將最弱邊的權(quán)值加到相鄰兩邊上，作為對去除最弱邊的補(bǔ)償。

圖1 鄰接三角形稀疏補(bǔ)償示意圖

分層稀疏化算法過程：

（1）對矩陣L中所有節(jié)點間的三角形，應(yīng)用稀疏化和懲罰過程。得到結(jié)果矩陣L。

（3）迭代過程的停止條件是矩陣L中節(jié)點之間的連接數(shù)減少到規(guī)定值。由算法2的迭代過程可得到不同稀疏程度的矩陣P的一個集合[P1P2P3P4P5P6... Pn]和最終的稀疏矩陣L～；由這兩者可用Krishnan等[8]中的算法1構(gòu)成一個分層預(yù)處理器F(x)，在求解線性方程Lx=b時，用F(x)代替Lx，減少方程的條件數(shù)以加快方程求解。

3 結(jié)果分析

3.1 對條件數(shù)的影響

對圖2應(yīng)用中的紋理剔除算法，在式中L求逆時應(yīng)用稀疏化處理。在矩陣迭代求逆過程中根據(jù)式（8）統(tǒng)計矩陣的條件數(shù)，圖3所示是處理前和處理后迭代過程中條件數(shù)的變化。由圖中信息可以看出，稀疏化處理可以減小條件數(shù)、加快迭代過程的收斂。

圖2 圖案樣本

圖3 分層稀疏化對矩陣條件數(shù)的影響

3.2 對處理時間的影響

收集500張如圖2所示圖片，從中提取出像素大小分別為64×64、256×256、1024×1024的圖案，組成測試處理時間的樣本庫。CPU為i5-4460，3.2GZ四核處理器，軟件運行環(huán)境是matlab2015a。以保邊濾波的結(jié)果求解過程為例測試分層稀疏化對圖像處理加速效果。即在式（3）中L求逆過程中，比較不應(yīng)用稀疏化處理和應(yīng)用稀疏化處理后迭代法解方程所需時間。

由表1中分層稀疏化前后平均求解時間的對比可以看出，對于像素行列數(shù)的圖案，分層稀疏化算法對求解時間能產(chǎn)生的影響有限，隨著圖案像素行列數(shù)增加，分層稀疏化算法對求解時間的影響越來越大。

表1 矩陣求逆平均時間比較Tab.1 Matrix inverse average time comparison

4 結(jié)束語

本文作者提出一種對大型線性方程的系數(shù)矩陣進(jìn)行預(yù)處理的方法。通過消去鄰域三角形中權(quán)重最小的邊、將權(quán)重補(bǔ)償?shù)狡渌麅蛇叢⒃谧酉到y(tǒng)重復(fù)迭代該過程的方式構(gòu)建分層預(yù)處理器。經(jīng)大量實例驗證，該算法對圖像處理中應(yīng)用廣泛的泊松方程求解有著很好的加速效果，可以大大減小系數(shù)矩陣的條件數(shù)，減少時間消耗。