双边滤波Matlab实现<The Bilateral Filter>

最新推荐文章于 2024-07-30 09:20:48 发布

沈子恒

最新推荐文章于 2024-07-30 09:20:48 发布

阅读量9.3k

点赞数 3

分类专栏：图像处理 | 目标跟踪 | OpenCV

图像处理 | 目标跟踪 | OpenCV 专栏收录该内容

68 篇文章 30 订阅

订阅专栏

原文：http://blog.csdn.net/lifeitengup/article/details/8902326#comments

双边滤波与一般的高斯滤波的不同就是：双边滤波既利用了位置信息<or 几何信息——高斯滤波只用了位置信息>又利用了像素信息来定义滤波窗口的权重。

像素值越接近，权重越大。双边滤波会去除图像的细节信息，又能保持边界。

对于彩色图像，像素值的接近与否不能使用RGB空间值，双边滤波的原始文献建议使用CIE颜色空间。

代码如下：

[plain]view plain copy 
   
 function resultI = BilateralFilt2(I,d,sigma)  
 %%%  
 %Author：LiFeiteng  
 %Version：1.0——灰色图像  Time：2013/05/01  
 %Version：1.1——灰色/彩色图像  Time：2013/05/02  2013/05/05  
 %d 半窗口宽度  
 I = double(I);  
 if size(I,3)==1  
     resultI = BilateralFiltGray(I,d,sigma);  
 elseif size(I,3)==3  
     resultI = BilateralFiltColor(I,d,sigma);  
 else   
     error('Incorrect image size')      
 end  
 end  
   
 function resultI = BilateralFiltGray(I,d,sigma)  
   
 [m n] = size(I);  
 newI = ReflectEdge(I,d);  
 resultI = zeros(m,n);  
 width = 2*d+1;  
 %Distance  
 D = fspecial('gaussian',[width,width],sigma(1));  
 S = zeros(width,width);%pix Similarity  
 h = waitbar(0,'Applying bilateral filter...');  
 set(h,'Name','Bilateral Filter Progress');  
 for i=1+d:m+d  
     for j=1+d:n+d  
         pixValue = newI(i-d:i+d,j-d:j+d);  
         subValue = pixValue-newI(i,j);  
         S = exp(-subValue.^2/(2*sigma(2)^2));  
         H = S.*D;  
         resultI(i-d,j-d) = sum(pixValue(:).*H(:))/sum(H(:));   
     end  
     waitbar(i/m);  
 end  
 close(h);  
 end  
   
 function resultI = BilateralFiltColor(I,d,sigma)  
 I = applycform(I,makecform('srgb2lab'));  
 [m n ~] = size(I);  
 newI = ReflectEdge(I,d);  
 resultI = zeros(m,n,3);  
 width = 2*d+1;  
 %Distance  
 D = fspecial('gaussian',[width,width],sigma(1));  
 % [X,Y] = meshgrid(-d:d,-d:d);  
 % D = exp(-(X.^2+Y.^2)/(2*sigma(1)^2));  
 S = zeros(width,width);%pix Similarity  
 h = waitbar(0,'Applying bilateral filter...');  
 set(h,'Name','Bilateral Filter Progress');  
 sigma_r = 100*sigma(2);  
 for i=1+d:m+d  
     for j=1+d:n+d  
         pixValue = newI(i-d:i+d,j-d:j+d,1:3);  
         %subValue = pixValue-repmat(newI(i,j,1:3),width,width);  
         dL = pixValue(:,:,1)-newI(i,j,1);  
         da = pixValue(:,:,2)-newI(i,j,2);  
         db = pixValue(:,:,3)-newI(i,j,3);  
         S = exp(-(dL.^2+da.^2+db.^2)/(2*sigma_r^2));  
         H = S.*D;  
         H = H./sum(H(:));  
         resultI(i-d,j-d,1) = sum(sum(pixValue(:,:,1).*H));   
         resultI(i-d,j-d,2) = sum(sum(pixValue(:,:,2).*H));      
         resultI(i-d,j-d,3) = sum(sum(pixValue(:,:,3).*H));      
     end  
     waitbar(i/m);  
 end  
 close(h);  
 resultI = applycform(resultI,makecform('lab2srgb'));  
 end  

其中newI = ReflectEdge(I,d); %对称地扩展边界，在原始图像I的边界处镜像映射像素值

[plain]view plain copy 
   
 function newI = ReflectEdge(I,d)  
 %Version：1.0——灰色图像  Time：2013/05/01  
 %Version：1.1——灰色/彩色图像  Time：2013/05/02  
 %考虑到实用性，决定不添加更多的边界处理选择，统一使用：reflect across edge  
   
 if size(I,3)==1  
     newI = ReflectEdgeGray(I,d);  
 elseif size(I,3)==3  
     newI = ReflectEdgeColor(I,d);  
 else   
     error('Incorrect image size')      
 end  
 end  
   
 function newI = ReflectEdgeGray(I,d)  
 [m n] = size(I);  
 newI = zeros(m+2*d,n+2*d);  
 %中间部分  
 newI(d+1:d+m,d+1:d+n) = I;  
 %上  
 newI(1:d,d+1:d+n) = I(d:-1:1,:);  
 %下  
 newI(end-d:end,d+1:d+n) = I(end:-1:end-d,:);  
 %左  
 newI(:,1:d) = newI(:,2*d:-1:d+1);  
 %右  
 newI(:,n+d+1:n+2*d) = newI(:,n+d:-1:n+1);  
 end  
   
 function newI = ReflectEdgeColor(I,d)  
 %扩展图像边界  
 [m n ~] = size(I);  
 newI = zeros(m+2*d,n+2*d,3);  
 %中间部分  
 newI(d+1:d+m,d+1:d+n,1:3) = I;  
 %上  
 newI(1:d,d+1:d+n,1:3) = I(d:-1:1,:,1:3);  
 %下  
 newI(end-d:end,d+1:d+n,1:3) = I(end:-1:end-d,:,1:3);  
 %左  
 newI(:,1:d,1:3) = newI(:,2*d:-1:d+1,1:3);  
 %右  
 newI(:,n+d+1:n+2*d,1:3) = newI(:,n+d:-1:n+1,1:3);  
 end  

测试用例：

[plain]view plain copy 
   
 img = imread('.\lena.tif');  
 %%img = imread('.\images\lena_gray.tif');  
 img = double(img)/255;  
 img = img+0.05*randn(size(img));  
 img(img<0) = 0; img(img>1) = 1;  
 %img = imnoise(img,'gaussian');  
 figure, imshow(img,[])  
 title('原始图像')  
 d = 6;  
 sigma = [3 0.1];  
 resultI = BilateralFilt2(double(img), d, sigma);  
   
 figure, imshow(resultI,[])  
 title('双边滤波后的图像')  

结果：

Reference：

1.C Tomasi, R Manduchi.Bilateral Filtering for Gray and Color Images, - Computer Vision, 1998.

像素值越接近，权重越大。双边滤波会去除图像的细节信息，又能保持边界。

对于彩色图像，像素值的接近与否不能使用RGB空间值，双边滤波的原始文献建议使用CIE颜色空间。

代码如下：

[plain]view plain copy 
   
 function resultI = BilateralFilt2(I,d,sigma)  
 %%%  
 %Author：LiFeiteng  
 %Version：1.0——灰色图像  Time：2013/05/01  
 %Version：1.1——灰色/彩色图像  Time：2013/05/02  2013/05/05  
 %d 半窗口宽度  
 I = double(I);  
 if size(I,3)==1  
     resultI = BilateralFiltGray(I,d,sigma);  
 elseif size(I,3)==3  
     resultI = BilateralFiltColor(I,d,sigma);  
 else   
     error('Incorrect image size')      
 end  
 end  
   
 function resultI = BilateralFiltGray(I,d,sigma)  
   
 [m n] = size(I);  
 newI = ReflectEdge(I,d);  
 resultI = zeros(m,n);  
 width = 2*d+1;  
 %Distance  
 D = fspecial('gaussian',[width,width],sigma(1));  
 S = zeros(width,width);%pix Similarity  
 h = waitbar(0,'Applying bilateral filter...');  
 set(h,'Name','Bilateral Filter Progress');  
 for i=1+d:m+d  
     for j=1+d:n+d  
         pixValue = newI(i-d:i+d,j-d:j+d);  
         subValue = pixValue-newI(i,j);  
         S = exp(-subValue.^2/(2*sigma(2)^2));  
         H = S.*D;  
         resultI(i-d,j-d) = sum(pixValue(:).*H(:))/sum(H(:));   
     end  
     waitbar(i/m);  
 end  
 close(h);  
 end  
   
 function resultI = BilateralFiltColor(I,d,sigma)  
 I = applycform(I,makecform('srgb2lab'));  
 [m n ~] = size(I);  
 newI = ReflectEdge(I,d);  
 resultI = zeros(m,n,3);  
 width = 2*d+1;  
 %Distance  
 D = fspecial('gaussian',[width,width],sigma(1));  
 % [X,Y] = meshgrid(-d:d,-d:d);  
 % D = exp(-(X.^2+Y.^2)/(2*sigma(1)^2));  
 S = zeros(width,width);%pix Similarity  
 h = waitbar(0,'Applying bilateral filter...');  
 set(h,'Name','Bilateral Filter Progress');  
 sigma_r = 100*sigma(2);  
 for i=1+d:m+d  
     for j=1+d:n+d  
         pixValue = newI(i-d:i+d,j-d:j+d,1:3);  
         %subValue = pixValue-repmat(newI(i,j,1:3),width,width);  
         dL = pixValue(:,:,1)-newI(i,j,1);  
         da = pixValue(:,:,2)-newI(i,j,2);  
         db = pixValue(:,:,3)-newI(i,j,3);  
         S = exp(-(dL.^2+da.^2+db.^2)/(2*sigma_r^2));  
         H = S.*D;  
         H = H./sum(H(:));  
         resultI(i-d,j-d,1) = sum(sum(pixValue(:,:,1).*H));   
         resultI(i-d,j-d,2) = sum(sum(pixValue(:,:,2).*H));      
         resultI(i-d,j-d,3) = sum(sum(pixValue(:,:,3).*H));      
     end  
     waitbar(i/m);  
 end  
 close(h);  
 resultI = applycform(resultI,makecform('lab2srgb'));  
 end  

其中newI = ReflectEdge(I,d); %对称地扩展边界，在原始图像I的边界处镜像映射像素值

[plain]view plain copy 
   
 function newI = ReflectEdge(I,d)  
 %Version：1.0——灰色图像  Time：2013/05/01  
 %Version：1.1——灰色/彩色图像  Time：2013/05/02  
 %考虑到实用性，决定不添加更多的边界处理选择，统一使用：reflect across edge  
   
 if size(I,3)==1  
     newI = ReflectEdgeGray(I,d);  
 elseif size(I,3)==3  
     newI = ReflectEdgeColor(I,d);  
 else   
     error('Incorrect image size')      
 end  
 end  
   
 function newI = ReflectEdgeGray(I,d)  
 [m n] = size(I);  
 newI = zeros(m+2*d,n+2*d);  
 %中间部分  
 newI(d+1:d+m,d+1:d+n) = I;  
 %上  
 newI(1:d,d+1:d+n) = I(d:-1:1,:);  
 %下  
 newI(end-d:end,d+1:d+n) = I(end:-1:end-d,:);  
 %左  
 newI(:,1:d) = newI(:,2*d:-1:d+1);  
 %右  
 newI(:,m+d+1:m+2*d) = newI(:,m+d:-1:m+1);  
 end  
   
 function newI = ReflectEdgeColor(I,d)  
 %扩展图像边界  
 [m n ~] = size(I);  
 newI = zeros(m+2*d,n+2*d,3);  
 %中间部分  
 newI(d+1:d+m,d+1:d+n,1:3) = I;  
 %上  
 newI(1:d,d+1:d+n,1:3) = I(d:-1:1,:,1:3);  
 %下  
 newI(end-d:end,d+1:d+n,1:3) = I(end:-1:end-d,:,1:3);  
 %左  
 newI(:,1:d,1:3) = newI(:,2*d:-1:d+1,1:3);  
 %右  
 newI(:,m+d+1:m+2*d,1:3) = newI(:,m+d:-1:m+1,1:3);  
 end  

测试用例：

[plain]view plain copy 
   
 img = imread('.\lena.tif');  
 %%img = imread('.\images\lena_gray.tif');  
 img = double(img)/255;  
 img = img+0.05*randn(size(img));  
 img(img<0) = 0; img(img>1) = 1;  
 %img = imnoise(img,'gaussian');  
 figure, imshow(img,[])  
 title('原始图像')  
 d = 6;  
 sigma = [3 0.1];  
 resultI = BilateralFilt2(double(img), d, sigma);  
   
 figure, imshow(resultI,[])  
 title('双边滤波后的图像')  