MatlabPython-WLS加权最小二乘滤波
最近看了論文Non-local Image Dehazing基于霧線的去霧算法。
其中算法的濾波用的是weighted least square(WLS)算法，論文全稱《Edge-Preserving Decompositions for Multi-Scale Tone and Detail Manipulation》，作者Z. Farbman等，發表在ACM SIGGRAPH 2007。
保邊濾波器可以是兩個矛盾的目標的結合體。對于一副輸入圖像g(N*M)，我們目標圖像u一方面我們希望其盡可能近似，與此同時除了在一些邊緣梯度變化比較大的地方外應該越平滑越好。數學模型為：
下標代表像素點空間位置，其中p表示像素的位置，ax和ay控制著不同位置上的平滑程度。
目標函數第一項代表輸入圖像和輸出圖像越相似越好，第二項是正則項，通過最小化的偏導，使得輸出圖像越平滑越好。平滑項的權重分別是和，依賴于輸入圖像，當輸入圖像的邊緣梯度變化比較大的時候，我們希望其約束小一些，以保留圖像的結構信息；當圖像的邊緣梯度變化很小的時候，這些細節信息我們認為其不重要，約束自然可以大一些，lambda是正則項參數，平衡兩者比重，越大圖像也就越平滑。?
上式寫成矩陣就是：
其中，Ax,Ay為以ax,ay為對角元素的對角矩陣，Dx，Dy為前向差分矩陣，和是后向差分算子，要使得（2）式去的最小值，u需滿足如下：
Lg為五點空間異構Laplacian矩陣。
論文中的去的平滑權重系數為：
輸入圖像亮度通道的log值（其實直接用原始圖像也是可行的），可以看出當的梯度比較大時，會隨著變小，否則反之。權重的這樣設計是很合理的，可以保留較大的邊緣，平滑不必要的細節。?其中?l?表示log，?ε一般取0.0001。
原作者在項目主頁中公布了源碼:
稀疏矩陣A就是公式(3)中的I+λLg。這段程序和前面的推導不同之處在于平滑權重先乘以了?λ，但是最中的結果是一致的：對角線上用1去減，其實就是加；對于其他四個對角線，我們推導的結果是帶有負號的平滑權重，然后乘以λ，其實就等價于直接乘以?λ。
wlsFilter.m
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矩陣的第一維度上做差分，也就是下面的行減去上面的行，得到(N-1)xM維的矩陣
dy = -lambda./(abs(dy).^alpha + smallNum);
dy = padarray(dy, [1 0], 'post');%在最后一行的后面補上一行0
dy = dy(:);%按列生成向量，就是Ay對角線上的元素構成的矩陣
dx = diff(L, 1, 2); %對L矩陣的第二維度做差分，也就是右邊的列減去左邊的列，得到Nx(M-1)的矩陣
dx = -lambda./(abs(dx).^alpha + smallNum);
dx = padarray(dx, [0 1], 'post');%在最后一列的后面添加一列0
dx = dx(:);%按列生成向量，對應上面Ax的對角線元素
% Construct a five-point spatially inhomogeneous Laplacian matrix
B(:,1) = dx;
B(:,2) = dy;
d = [-r,-1];
A = spdiags(B,d,k,k);%把dx放在-r對應的對角線上，把dy放在-1對應的對角線上
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只有五個對角線上有非0元素
% Solve
OUT = A\IN(:);
OUT = reshape(OUT, r, c);
?
demo_wlsFilter.m
img = imread('C:\Users\x\Desktop\15.jpg');
IN = double(img);
OUT = IN;
%OUT = wlsFilter(IN, 0.5);
if length(size(img))>2
? ? for i=1:3
? ? ? ? IN(:,:,i) = IN(:,:,i)/max(IN(:));
? ? ? ? OUT(:,:,i) = wlsFilter(IN(:,:,i),0.35);
? ? end
end
figure('Position',[50,50, size(IN,2)*2 , size(IN,1)]);
subplot(1,3,1); imshow(IN);title('input')
subplot(1,3,2); imshow(OUT); title('output')
imwrite(OUT,['C:\Users\x\Desktop\', int2str(16),'.jpg']);
平滑結果：
Python實現
由于numpy的矩陣很容易溢出出現MemoryError，而A的大小很容易超過10萬行&列，所以圖片必須壓縮。測試代碼用的是（128x128）大小的圖片。
還有一點就是python的spdiags與matlab的參數是不一樣的。spdiags(B, diags, k, k)這里的B要用matlab的B參數的轉置。另外python沒有\除，所以矩陣a必須先求逆矩陣，由于矩陣是16384x16384大小，所以代碼要跑3分鐘。。。
import numpy as np import cv2 #import math from scipy.sparse import spdiagsdef wlsFilter(img, Lambda, alpha=1.2, eps=0.0001):L = np.log(img+0.0001)row, cols = img.shape[:2]k = row * cols#對L矩陣的第一維度上做差分，也就是下面的行減去上面的行，得到(N-1)xM維的矩陣dy = np.diff(L, 1, 0)dy = -Lambda/(np.power(np.abs(dy),alpha)+eps)#在最后一行的后面補上一行0dy = np.pad(dy, ((0,1),(0,0)), 'constant')#按列生成向量，對應上面Ay的對角線元素dy = dy.Tdy = dy.reshape(-1,1)#對L矩陣的第二維度上做差分，也就是右邊的列減去左邊的列，得到Nx(M-1)的矩陣dx = np.diff(L, 1, 1)dx = -Lambda/(np.power(np.abs(dx),alpha)+eps)#在最后一列的后面補上一行0dx = np.pad(dx, ((0,0),(0,1)), 'constant')#按列生成向量，對應上面Ay的對角線元素dx = dx.Tdx = dy.reshape(-1,1)#構造五點空間非齊次拉普拉斯矩陣B = np.hstack((dx,dy))B = B.Tdiags = np.array([-row, -1])#把dx放在-row對應的對角線上，把dy放在-1對應的對角線上A = spdiags(B, diags, k, k).toarray()e = dxw = np.pad(dx, ((row,0),(0,0)), 'constant')w = w[0:-row]s = dyn = np.pad(dy, ((1,0),(0,0)), 'constant')n = n[0:-1]D = 1-(e+w+s+n)D = D.T#A只有五個對角線上有非0元素diags1 = np.array([0])A = A + np.array(A).T + spdiags(D, diags1, k, k).toarray()im = np.array(img)p,q = im.shape[:2]g = p*qim = np.reshape(im,(g,1))a = np.linalg.inv(A)out = np.dot(a,im)out = np.reshape(out,(row, cols))return outimg = cv2.imread(r'C:\Users\x\Desktop\k2.jpg', cv2.IMREAD_ANYCOLOR)m = np.double(img)b, g, r = cv2.split(m)ib = np.array(b) p1,q1 = ib.shape[:2] h1 = p1*q1 ib = np.reshape(ib,(h1,1)) b = b/np.max(ib)ig = np.array(g) p2,q2 = ig.shape[:2] h2 = p2*q2 ig = np.reshape(ig,(h2,1)) g = g/np.max(ig)ir = np.array(r) p3,q3 = ir.shape[:2] h3 = p3*q3 ir = np.reshape(ir,(h3,1)) r = r/np.max(ir)wls1 = wlsFilter(b,1) wls2 = wlsFilter(g,1) wls3 = wlsFilter(r,1) wls = cv2.merge([wls1, wls2, wls3])cv2.imshow('image', img) cv2.imshow('filter', wls) cv2.imwrite(r'C:\Users\x\Desktop\18.jpg', wls) cv2.waitKey(0) cv2.destroyAllWindows()?
總結
以上是生活随笔為你收集整理的MatlabPython-WLS加权最小二乘滤波的全部內容，希望文章能夠幫你解決所遇到的問題。
