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