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

墨蓝 2022-08-08 15:52 17阅读 0赞

## 引言 ##

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

## 算法 ##

设计一个保边滤波器可以看做是两个矛盾的目标的结合体。对于一副输入图像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][forward difference operators]）。那就是我们的  Lg 是一个五点空间异性拉普拉斯矩阵（ [five-point spatially inhomogeneous Laplacian  
matrix][five-point spatially inhomogeneous Laplacian matrix]）。简单解释一下异性，异性因为  Lg 中  Ax 是个变量，随空间位置改变。若它是一个常数，即随空间位置不变，那么这样的矩阵就称为五点空间同性拉普拉斯矩阵。更多内容可以参考 [《Efficient Preconditioning of Laplacian Matrices for Computer Graphics》][Efficient Preconditioning of Laplacian Matrices for Computer Graphics]一文，里面介绍了比较多的拉普拉斯矩阵的数学基础和应用。  
这个算法的核心还是生成拉普拉斯这个矩阵。  
先准备一下数据，一会儿再来填充一下这个矩阵。  
第一步根据图像  l 的梯度值计算相邻像素之间的平滑权重：

% 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一会儿将被填充到非主对角线位置，这些元素都是为负的。  
接下来就是构造拉普拉斯矩阵了，这里的拉普拉斯是一个对称，只有少数几个对角线有元素，其余为零的稀疏矩阵。

% 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是上（北）下（南）左（西）右（东）四个方位。  
看最终生成的一副拉普拉斯矩阵图吧。  
![这里写图片描述][20150921204802671]  
这张图中可以看出每一行元素之和都为0。其中紧靠主对角线元素的两个对角线填充的是dy元素,比较远的对角线填充的是dx元素，这样拉普拉斯矩阵处理的就是二维图像了。  
完整的代码如下：

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);%转换为矩阵

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

## 关于拉普拉斯矩阵 ##

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

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).][Link 1]  
[Krishnan D, Fattal R, Szeliski R. 2013. Efficient preconditioning of laplacian matrices for computer graphics\[J\]. ACM Transactions on Graphics (Proc. SIGGRAPH), 32,4][Efficient Preconditioning of Laplacian Matrices for Computer Graphics]

## 转载请保留以下信息 ##

<table> 
 <thead> 
  <tr> 
   <th>作者</th> 
   <th align="center">日期</th> 
   <th align="right">联系方式</th> 
  </tr> 
 </thead> 
 <tbody> 
  <tr> 
   <td>风吹夏天</td> 
   <td align="center">2015年9月21日</td> 
   <td align="right">wincoder@qq.com</td> 
  </tr> 
 </tbody> 
</table>

[20150919161836329]: /images/20220731/917e48e2e2f942fc884d49aee2fc65f0.png
[Link 1]: http://www.cs.huji.ac.il/~danix/epd/
[forward difference operators]: https://en.wikipedia.org/wiki/Finite_difference
[five-point spatially inhomogeneous Laplacian matrix]: https://en.wikipedia.org/wiki/Laplacian_matrix
[Efficient Preconditioning of Laplacian Matrices for Computer Graphics]: http://www.cs.huji.ac.il/~raananf/projects/hsc/
[20150921204802671]: /images/20220731/d6df74f0a78d48559d2198f5ad3c882f.png
[Link 2]: http://www.cs.huji.ac.il/~yweiss/Colorization/
[Link 3]: https://www.cs.jhu.edu/~misha/Fall07/Papers/Perez03.pdf
[HDR]: http://www.cs.huji.ac.il/~danix/hdr/