引言
陆陆续续在计算摄影学接触了不少保边滤波器，其重要性自不必说，可以用在图像的增强，图像抽象画，高动态范围图像压缩，图像色调映射等。
今天介绍的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 = AIN(:); %公式（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