clear all;close all;

%% 读取并显示原始图像
I=imread('lena.bmp');
figure;imshow(I);

%% 加噪音并显示
I=imnoise(I,'gaussian',0,0.001);
figure;imshow(I);

%% 高斯滤波
[m,n]=size(I);
winsize=3;                      % wisize----模板大小必须是奇数
padwidth=(winsize-1)/2;         % padwidth----模板宽度
sigma=0.5;                      % 高斯模板----体现空间距离相关性
gaus_filter = fspecial('gaussian',winsize,winsize/6);

for i=2:m-1
    for j=2:n-1
        %weight:权重模板-----体现灰度相关性
        temp=I(i-padwidth:i+padwidth,j-padwidth:j+padwidth);
        temp=double(temp);
        a=temp.*temp;
        weight=exp(-abs(a-a(padwidth+1,padwidth+1))/2/sigma/sigma);
        weight=-(temp-temp(padwidth+1,padwidth+1))/255+1;
        coef=sum(sum(weight));
        weight=weight/coef;
        
        bilateral_filter=gaus_filter.*weight;
       
       New=temp.*bilateral_filter;
       Newpix=sum(sum(New));
       J(i,j)=Newpix;
    end
end
J=uint8(J);
J=histeq(J);  %亮度偏暗，需均衡化
figure;
imshow(J);