加權最小二乘濾波WLS（weighted least squares）加上雙邊濾波，引導濾波是三種較為經典的邊緣保持性濾波演算法，該演算法最早見於論文：《Edge-Preserving Decompositions for Multi-Scale Tone and Detail Manipulation》中，原作者專案主頁：http://www.cs.huji.ac.il/~danix/epd/，本篇進行總結和測試。

加權最下二乘濾波原理：

作者提出該演算法的初衷是，基於雙邊濾波的演算法無法在多尺度上提取到很好的細節資訊，並可能出現偽影。加權最小二乘濾波目的即是使得結果影象u 與原始影象 p經過平滑後儘量相似，但是在邊緣部分儘量保持原狀，用數學表達出來即為：

  （1）

其中，ax,ay為權重係數。目標函式第一項(up−gp)2代表輸入影象u和輸出影象g越相似越好；第二項是正則項，通過最小化u的偏導，使得輸出影象g越平滑越好。

上式可以改寫為矩陣形式：

                 （2）

其中，Ax,Ay為以ax,ay為對角元素的對角矩陣，Dx，Dy為前向差分矩陣，DTx和DTy要使得（2）式去的最小值，u需滿足如下：

                                                                     （3）

其中，，作者去的平滑權重係數為：

其中 l 表示log， ε一般取0.0001。式（3）求出μ為：

，                                                     （4）

當所選區域連續是，平滑權重係數可以近似為：ax ≈ ay ≈ a,那麼：

                                                                     （5）

加權最小二乘濾波實現：

原作者在專案主頁中公佈了原始碼，這裡po一下，便於理解：

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);  %對L矩陣的第一維度上做差分
dy = -lambda./(abs(dy).^alpha + smallNum);
dy = padarray(dy, [1 0], 'post'); %在最後一行後面補一行0
dy = dy(:); %按列生成向量，就是Ay對角線上元素構成的矩陣，下同

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


% 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);

% Solve
OUT = A\IN(:);
OUT = reshape(OUT, r, c);
 

這裡A+A′構造的是拉普拉斯非主對角線元素，D是主對角線元素。n,s,w,e是上（北）下（南）左（西）右（東）四個方位。 最終生成的一副拉普拉斯矩陣圖：

圖中每一行元素之和都為0。其中緊靠主對角線元素的兩個對角線填充的是dy元素,比較遠的對角線填充的是dx元素，這樣拉普拉斯矩陣處理的就是二維影象了。 

%% 加權最小二乘濾波測試函式
clc,close all,clear all;
RGB = imread('flower.png'); 
% if length(size(RGB))>2
%     I = double(rgb2gray(RGB));
% else
%     I=double(RGB);
% end
I=double(RGB);
res=I;
if length(size(RGB))>2
    for i=1:3
        I(:,:,i) = I(:,:,i)./max(I(:));
        res(:,:,i) = wlsFilter(I(:,:,i));
    end
end
figure,
subplot(211),imshow(I),title('原圖');
subplot(212), imshow(res),title('wls-output');

參考：

http://blog.csdn.net/bluecol/article/details/48576253