基于加权最小二乘法的保边缘平滑滤波器（WLS）

引言

陆陆续续在计算摄影学接触了不少保边滤波器，其重要性自不必说，可以用在图像的增强，图像抽象画，高动态范围图像压缩，图像色调映射等。 
 
今天介绍的WLS(最小二乘滤波器）即使其中一种，论文全称《Edge-Preserving Decompositions for Multi-Scale Tone and Detail Manipulation》，作者Z. Farbman等，发表在ACM SIGGRAPH 2007。 
这篇文章是基于最优化理论，代码是公开，点此下载。当然有更好的滤波器不断出现。但是这篇文章核心代码只有十几行左右，数学功底深厚，方法比较去值得研究和学习之用。这篇文章只是我看此论文一个笔录，不正确之处欢迎指正。

算法

设计一个保边滤波器可以看做是两个矛盾的目标的结合体。对于一副输入图像 g ，我们目标图像 u 一方面我们希望其尽可能近似 g ，与此同时 u 除了在 g 一些边缘梯度变化比较大的地方外应该越平滑越好。形式上，我们寻求最小化下列目标函数的解。 

∑p((up−gp)2+λ(ax,p(g)(∂u∂x)2)+ay,p(g)(∂u∂y)2))(1)

这里，下标 p 代表像素点空间位置。目标函数第一项 (up−gp)2 代表输入图像 u 和输出图像 g 越相似越好，第二项是正则项，通过最小化 u 的偏导，使得输出图像 g 越平滑越好。平滑项的权重分别是 ax 和 ay ，依赖于输入图像 g ，想想也是这样的嘛，当输入图像的边缘梯度变化比较大的时候，我们希望其约束小一些，以保留图像的结构信息，当图像的边缘梯度变化很小的时候，这些细节信息我们认为其不重要，约束自然可以大一些， λ 是正则项参数，平衡两者比重， λ 越大图像也就越平滑。 
将上式写成矩阵形式（因为矩阵比较清晰简洁，很好用的工具，^_^）： 

(u−g)T(u−g)+λ(uTDTxAxDxu+uTDTyAyDyu)(2)

这里， Ax 和 Ay 是包含平滑权重 ax(g),ay(g) 的对角矩阵，一会儿我们将看到其是输入图像 g 梯度的函数。 Dx和Dy 是离散差分算子。 
这里 u 和 g 都会转换成向量形式，对上式求导（比较简单的一步）可得到下列线性方程组， 

(I+λLg)u=g(3)

这里 Lg=DTxAxDx+DTyAyDy 。 
还有平滑项权重的设置，作者如下定义 

ax,p(g)=(|∂l∂x(p)|α+ε)−1(4)

ay,p(g)=(|∂l∂y(p)|α+ε)−1(5)

l  是输入图像亮度通道的log值（其实直接用原始图像也是可行的）,可以看出当 l 的梯度比较大时， ax,p(g),ay,p(g) 会随着变小，否则反之。权重的这样设计是很合理的，可以保留较大的边缘，平滑不必要的细节。 
注意到 Dx 和 Dy 是前向算子（ forward difference operators ）， DTx 和 DTy 是后向差分算子（ backward difference opeattors ）。那就是我们的 Lg 是一个五点空间异性拉普拉斯矩阵（ five-point spatially inhomogeneous Laplacian 
matrix ）。简单解释一下异性，异性因为 Lg 中 Ax 是个变量，随空间位置改变。若它是一个常数，即随空间位置不变，那么这样的矩阵就称为五点空间同性拉普拉斯矩阵。更多内容可以参考 《Efficient Preconditioning of Laplacian Matrices for Computer Graphics》 一文，里面介绍了比较多的拉普拉斯矩阵的数学基础和应用。 
这个算法的核心还是生成拉普拉斯这个矩阵。 
先准备一下数据，一会儿再来填充一下这个矩阵。 
第一步根据图像 l 的梯度值计算相邻像素之间的平滑权重：

<pre name="code" class="plain">% Compute affinities between adjacent pixels based on gradients of L
dy = diff(L, 1, 1);
dy = -lambda./(abs(dy).^alpha + smallNum);
dy = padarray(dy, [1 0], 'post');
dy = dy(:);

dx = diff(L, 1, 2); 
dx = -lambda./(abs(dx).^alpha + smallNum);
dx = padarray(dx, [0 1], 'post');
dx = dx(:);

这一步基本上就是按照公式（4）（5）编写的，只是lambda前面多了一个负号， dx 和 dy 一会儿将被填充到非主对角线位置，这些元素都是为负的。 
接下来就是构造拉普拉斯矩阵了，这里的拉普拉斯是一个对称，只有少数几个对角线有元素，其余为零的稀疏矩阵。

<pre name="code" class="plain">% Construct a five-point spatially inhomogeneous Laplacian matrix
B(:,1) = dx;
B(:,2) = dy;
d = [-r,-1];
A = spdiags(B,d,k,k);

e = dx;
w = padarray(dx, r, 'pre'); w = w(1:end-r);
s = dy;
n = padarray(dy, 1, 'pre'); n = n(1:end-1);

D = 1-(e+w+s+n);
A = A + A' + spdiags(D, 0, k, k);

这里 A+A′ 构造的是拉普拉斯的非主对角线元素， D 是主对角线元素。 n,s,w,e 是上（北）下（南）左（西）右（东）四个方位。 
看最终生成的一副拉普拉斯矩阵图吧。 
 
这张图中可以看出每一行元素之和都为0。其中紧靠主对角线元素的两个对角线填充的是 dy 元素,比较远的对角线填充的是 dx 元素，这样拉普拉斯矩阵处理的就是二维图像了。 
完整的代码如下：

<pre name="code" class="plain">function OUT = wlsFilter(IN, lambda, alpha, L)
%WLSFILTER Edge-preserving smoothing based on the weighted least squares(WLS) 
%   optimization framework, as described in Farbman, Fattal, Lischinski, and
%   Szeliski, "Edge-Preserving Decompositions for Multi-Scale Tone and Detail
%   Manipulation", ACM Transactions on Graphics, 27(3), August 2008.
%
%   Given an input image IN, we seek a new image OUT, which, on the one hand,
%   is as close as possible to IN, and, at the same time, is as smooth as
%   possible everywhere, except across significant gradients in L.
%
%
%   Input arguments:
%   ----------------
%     IN              Input image (2-D, double, N-by-M matrix). 
%       
%     lambda          Balances between the data term and the smoothness
%                     term. Increasing lambda will produce smoother images.
%                     Default value is 1.0
%       
%     alpha           Gives a degree of control over the affinities by non-
%                     lineary scaling the gradients. Increasing alpha will
%                     result in sharper preserved edges. Default value: 1.2
%       
%     L               Source image for the affinity matrix. Same dimensions
%                     as the input image IN. Default: log(IN)
% 
%
%   Example 
%   -------
%     RGB = imread('peppers.png'); 
%     I = double(rgb2gray(RGB));
%     I = I./max(I(:));
%     res = wlsFilter(I, 0.5);
%     figure, imshow(I), figure, imshow(res)
%     res = wlsFilter(I, 2, 2);
%     figure, imshow(res)

if(~exist('L', 'var')),
    L = log(IN+eps);
end

if(~exist('alpha', 'var')),
    alpha = 1.2;
end

if(~exist('lambda', 'var')),
    lambda = 1;
end

smallNum = 0.0001;

[r,c] = size(IN);
k = r*c;

% Compute affinities between adjacent pixels based on gradients of L
dy = diff(L, 1, 1);
dy = -lambda./(abs(dy).^alpha + smallNum);
dy = padarray(dy, [1 0], 'post');
dy = dy(:); %公式（5）

dx = diff(L, 1, 2); 
dx = -lambda./(abs(dx).^alpha + smallNum);
dx = padarray(dx, [0 1], 'post');
dx = dx(:); %公式（4）


% Construct a five-point spatially inhomogeneous Laplacian matrix
B(:,1) = dx;
B(:,2) = dy;
d = [-r,-1];
A = spdiags(B,d,k,k);

%构造主对角线
e = dx;
w = padarray(dx, r, 'pre'); w = w(1:end-r);
s = dy;
n = padarray(dy, 1, 'pre'); n = n(1:end-1);

D = 1-(e+w+s+n); %再加上单位矩阵，这里元素都为负数，先取反
A = A + A' + spdiags(D, 0, k, k); %A+A'为非主对角线元素

% Solve
OUT = A\IN(:); %公式（3）
OUT = reshape(OUT, r, c);%转换为矩阵

算法介绍完毕，关于其应用可以参考论文的主页，见末尾。

关于拉普拉斯矩阵

拉普拉斯矩阵在计算机图形和计算摄影学经常会遇到的，比如说灰度图像着色问题，泊松图像编辑，HDR，还有保边滤波器的设计，当然在计算机图形学中的也有诸多应用，所以掌握其解法还是非常有必要的。 
它经常会出现在如下的最优化问题当中 

F(x)=∑i∈I[ui(xi−yi)2+∑j∈Niwij(xi−xj−zij)2]

第一项是数据项，度量 x 与输入向量 y 的差异，第二项是平滑项，度量每一个变量 xi 与其所在局部窗口中的邻近像素 xj,j∈Ni 的差异，相对于 Zij 。一般情况下这个局部窗口我们取当前像素点的四邻域或者八邻域。 ui 和 wij 都是一些权重信息，与所研究的问题相关，都是非负的。当 ui 和 wij 是常数时候，这个问题就是空间同性拉普拉斯，否则就是空间异性拉普拉斯问题了。

参考资料

FARBMAN, Z., FATTAL, R., LISCHINSKI, D., AND SZELISKI, R. 2008. Edge-preserving decompositions for multi-scale tone and detail manipulation. ACM Transactions on Graphics (Proc. SIGGRAPH) 27, 3 (August). 
Krishnan D, Fattal R, Szeliski R. 2013. Efficient preconditioning of laplacian matrices for computer graphics[J]. ACM Transactions on Graphics (Proc. SIGGRAPH), 32,4

转载请保留以下信息

作者	日期	联系方式
风吹夏天	2015年9月21日	wincoder@qq.com