局部边缘保持滤波(LEP)高动态范围图像HDR压缩 matlab程序（一）

本文为转载，原博客地址为：http://blog.csdn.net/majinlei121/article/details/50420927

本文程序根据论文"Gu B, Li W, Zhu M, et al. Local edge-preserving multiscale decomposition for high dynamic range image tone mapping[J]. Image Processing, IEEE Transactions on, 2013, 22(1): 70-79."编写，主要为论文中提出的一种LEP局部边缘保持滤波程序，高动态图像压缩的程序在下一篇博客里。

下面为LEP函数程序：

[cpp] view plain copy

1. function out = LocalWls_HDR(I, alpha, beta, r)  
2.   
3. if ~exist('alpha','var')  
4.     alpha = 0.1;  
5. end  
6. if ~exist('beta','var')  
7.     beta = 1;  
8. end  
9. if ~exist('r','var')  
10.     r = 4;  
11. end  
12.   
13. [hei, wid] = size(I);  
14. N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.  
15.   
16. mean_I = boxfilter(I, r) ./ N;  
17. mean_II = boxfilter(I.*I, r) ./ N;  
18.   
19. var_I = mean_II - mean_I .* mean_I;  
20.   
21. I_x = diff(I,1,2);  
22. I_x = padarray(I_x, [0 1 0], 'post');%进行列方向差分，求dx  
23. I_y = diff(I,1,1);  
24. I_y = padarray(I_y, [1 0 0], 'post');%进行行方向差分，求dy  
25. I_gradient=abs(I_x+I_y);  
26. I_gradient=I_gradient.^(2-beta);  
27. I_gradient=alpha*boxfilter(I_gradient,r)./N;  
28.   
29. a = var_I ./ (var_I + I_gradient+0.00000001);   
30. b = mean_I - a .* mean_I;   
31.   
32. mean_a = boxfilter(a, r) ./ N;  
33. mean_b = boxfilter(b, r) ./ N;  
34.   
35. out = mean_a .* I + mean_b;   
36. end  

其中用到的boxfilter()函数来自引导滤波，为：

[cpp] view plain copy

1. function imDst = boxfilter(imSrc, r)  
2.   
3. %   BOXFILTER   O(1) time box filtering using cumulative sum  
4. %  
5. %   - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));  
6. %   - Running time independent of r;   
7. %   - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);  
8. %   - But much faster.  
9.   
10. [hei, wid] = size(imSrc);  
11. imDst = zeros(size(imSrc));  
12.   
13. %cumulative sum over Y axis  
14. imCum = cumsum(imSrc, 1);  
15. %difference over Y axis  
16. imDst(1:r+1, :) = imCum(1+r:2*r+1, :);  
17. imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);  
18. imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);  
19.   
20. %cumulative sum over X axis  
21. imCum = cumsum(imDst, 2);  
22. %difference over X axis  
23. imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);  
24. imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);  
25. imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);  
26. end  

运行程序为：

[cpp] view plain copy

1. clear all  
2. I = double(imread('..\images\stripes.png')) / 255;  
3.   
4. alpha=0.1;  
5. beta=1;  
6. r=2;  
7.   
8. out = LocalWls_HDR(I, alpha, beta, r);  
9.   
10. figure();imshow([I, out], [0, 1]);title('input,LEP');  
11. row=155;  
12. figure;plot(1:size(I,2), I(row,:,1),1:size(I,2), out(row,:,1), 'LineWidth', 2);  
13. legend('input','LEP');  

效果图为：