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

田俊杰，趙祖燁，張軍飛

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

0 引言

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

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

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

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

1 算法應用場景

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

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

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

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

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

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

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

2 算法分析

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

由于L的半正定型，能量值總是非負數(shù)。矩陣L的特征向量對應的能量值是特征向量的對應的特征值，由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和細節(jié)點F，通過迭代過程分層稀疏化矩陣。矩陣L可以劃分為細節(jié)點C之間連接LCC，粗節(jié)點F之間連接LFF、細節(jié)點C與粗節(jié)點F之間連接LFC：

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

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

圖1 鄰接三角形稀疏補償示意圖

分層稀疏化算法過程：

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

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

3 結(jié)果分析

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

對圖2應用中的紋理剔除算法，在式中L求逆時應用稀疏化處理。在矩陣迭代求逆過程中根據(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求逆過程中，比較不應用稀疏化處理和應用稀疏化處理后迭代法解方程所需時間。

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

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

4 結(jié)束語

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