9月12日【边缘检测、锐化】

利用拉普拉斯算子对模糊图像进行增强

clear all;   %利用拉普拉斯算子对模糊图像进行增强
I=imread('rice.png');
subplot(1,2,1);imshow(I);
title('原始图像');
I=double(I);  %转换数据类型为double双精度型
H=[0 1 0,1 -4 1,0 1 0];  %拉普拉斯算子
J=conv2(I,H,'same');  %用拉普拉斯算子对图像进行二维卷积运算
K=I-J; 
subplot(1,2,2),imshow(K,[])
title('锐化滤波处理')


对图像进行梯度法锐化
clear all;
[I,map]=imread('trees.tif');
subplot(2,2,1);imshow(I);
title('原始图像');
I=double(I);                    %数据类型转换
[IX,IY]=gradient(I);            %梯度
gm=sqrt(IX.*IX+IY.*IY);
out1=gm;
subplot(2,2,2);imshow(out1,map);
title('梯度值');
out2=I;
J=find(gm>=15);      %阈值处理
out2(J)=gm(J);
subplot(2,2,3);imshow(out2,map);
title('对梯度值加阈值');
out3=I;
J=find(gm>=20);      %阈值黑白化
out3(J)=255;         %设置为白色
K=find(gm<20);       %阈值黑白化
out3(K)=0;           %设置为黑色
subplot(2,2,4);imshow(out3,map);  
title('二值化处理')

分别使用sobel、prewitt、log、canny算子对图像进行锐化

I=imread('coins.png');%读入图像
subplot(3,2,1),
imshow(I);
title('原始图像');%显示原图像
H=fspecial('sobel'); %应用sobel算子锐化图像
I2=filter2(H,I); %sobel算子滤波锐化
subplot(3,2,2);
imshow(I2); %显示sobel算子锐化图像
title('sobel算子锐化图像');
H1=fspecial('prewitt');%应用prewitt算子锐化图像
I3=filter2(H1,I);%prewitt算子滤波锐化
subplot(3,2,3);imshow(I3); %显示prewitt算子锐化图像
title('prewitt算子锐化图像');
H2=fspecial('log'); %应用log算子锐化图像
I4=filter2(H2,I); %log算子滤波锐化
subplot(3,2,4);
imshow(I4);%显示log算子锐化图像
title('log算子锐化图像');
H3=edge(I,'canny',[0.04 0.10],1.5); %应用log算子锐化图像
% I5=filter2(H3,I); %log算子滤波锐化
subplot(3,2,5);
imshow(H3,[]);%显示canny算子锐化图像
title('canny算子锐化图像');

低通滤波器

图像在传递过程中，由于噪声主要集中在高频部分，为去除噪声改善图像质量，滤波器采用低通滤波器H（u，v）来抑制高频成分，通过低频成分，然后再进行傅里叶逆变换获得滤波图像，就可达到平滑图像的目的。由卷积定理，低通滤波器数学表达式为：
G(u,v)=F(u,v)*H(u,v),其中，F(u,v)为含噪声的原图像的傅里叶变换域，H（u,v）为传递函数，G(u,v)为经过低通滤波后输出图像的傅里叶变换。
1，理想的低通滤波器。
I = imread('tire.tif');
[f1,f2] = freqspace(size(I),'meshgrid'); %生成频率序列矩阵
Hd = ones(size(I));   
r = sqrt(f1.^2 + f2.^2); 
Hd(r>0.1) = 0;   %构造滤波器
Y=fft2(double(I));   
Y=fftshift(Y);   
Ya=Y.*Hd;   %滤波
Ya=ifftshift(Ya);   
Ia01=ifft2(Ya);
Hd(r>0.2) = 0;   %构造滤波器
Y=fft2(double(I));  
Y=fftshift(Y);   
Ya=Y.*Hd;   
Ya=ifftshift(Ya);  
Ia02=ifft2(Ya);
Hd(r>0.5) = 0;   %构造滤波器
Y=fft2(double(I));  
Y=fftshift(Y); 
Ya=Y.*Hd;   
Ya=ifftshift(Ya);  
Ia05=ifft2(Ya);
subplot(2,2,1),
imshow(I),
title('原始图像')
subplot(2,2,2),
imshow(uint8(Ia01)),
title('r为0.1时')
subplot(2,2,3),
imshow(uint8(Ia02)),
title('r为0.2时')
subplot(2,2,4),
imshow(uint8(Ia05)),
title('r为0.5时')


2.巴特沃斯低通滤波器。

I = imread('cell.tif');
[f1,f2] =freqspace(size(I),'meshgrid');
D=0.4;     %截止频率
n=1;
Hd = ones(size(I)); 
r = sqrt(f1.^2 + f2.^2);
for i=1:size(I,1)
    for j=1:size(I,2)
        t=r(i,j)/(D*D);
        Hd(i,j) = 1/(t^n+1);%构造滤波函数
    end
end
B=fft2(double(I)); 
B=fftshift(B);
Ba=B.*Hd;
Ba=ifftshift(Ba);
Ia1=ifft2(Ba);
n=2;    
Hd = ones(size(I)); 
r = sqrt(f1.^2 + f2.^2);
for i=1:size(I,1)
    for j=1:size(I,2)
        t=r(i,j)/(D*D);
        Hd(i,j) = 1/(t^n+1);%构造滤波函数
    end
end
B=fft2(double(I)); 
B=fftshift(B);%
Ba=B.*Hd;%
Ba=ifftshift(Ba);%
Ia2=ifft2(Ba);%
n=6;    
Hd = ones(size(I)); 
r = sqrt(f1.^2 + f2.^2);
for i=1:size(I,1)
    for j=1:size(I,2)
        t=r(i,j)/(D*D);
        Hd(i,j) = 1/(t^n+1);%构造滤波函数
    end
end
B=fft2(double(I)); 
B=fftshift(B);
Ba=B.*Hd;
Ba=ifftshift(Ba);
Ia6=ifft2(Ba);
subplot(2,2,1),
imshow(I),
title('原始图像');
subplot(2,2,2),
imshow(uint8(Ia1)),
title('n为10时的滤波效果');
subplot(2,2,3),
imshow(uint8(Ia2)),
title('n为13时的滤波效果');
subplot(2,2,4),
imshow(uint8(Ia6)),
title('n为18时的滤波效果');


3.梯形低通滤波器。

4指数低通滤波器。
I = imread('gantrycrane.png');
[f1,f2] = freqspace(size(I),'meshgrid');
D=10/size(I,1); %D为10时
Hd = ones(size(I)); 
r = f1.^2 + f2.^2;
for i=1:size(I,1)
    for j=1:size(I,2)
        t=r(i,j)/(D*D);
        Hd(i,j) = exp(-t);
    end
end
E=fft2(double(I));
E=fftshift(E);
Ea=E.*Hd;
Ea=ifftshift(Ea);
Ia10=ifft2(Ea);
D=40/size(I,1); %D为40时
Hd = ones(size(I)); 
r = f1.^2 + f2.^2;
for i=1:size(I,1)
    for j=1:size(I,2)
        t=r(i,j)/(D*D);
        Hd(i,j) = exp(-t);
    end
end
E=fft2(double(I));
E=fftshift(E);
Ea=E.*Hd;
Ea=ifftshift(Ea);
Ia40=ifft2(Ea);
D=100/size(I,1); %D为100时
Hd = ones(size(I)); 
r = f1.^2 + f2.^2;
for i=1:size(I,1)
    for j=1:size(I,2)
        t=r(i,j)/(D*D);
        Hd(i,j) = exp(-t);
    end
end
E=fft2(double(I));
E=fftshift(E);
Ea=E.*Hd;
Ea=ifftshift(Ea);
Ia100=ifft2(Ea);
subplot(2,2,1),
imshow(I),
title('原始图像') %图像显示
subplot(2,2,2),
imshow(uint8(Ia10)),
title('D为10时');
subplot(2,2,3),
imshow(uint8(Ia40)),
title('D为40时');
subplot(2,2,4),
imshow(uint8(Ia100)),
title('D为100时');


各种低通滤波器对图像滤波

clear all;

[I,map]=imread('coins.png');     %读取图像

noisy=imnoise(I,'gaussian',0.02);

[M,N]=size(I);

F=fft2(noisy);

fftshift(F);

Dcut=100;

D0=150;

D1=250;

for u=1:M

    for v=1:N

        D(u,v)=sqrt(u^2+v^2);

        BUTTERH(u,v)=1/(1+(sqrt(2)-1)*(D(u,v)/Dcut)^2);

        EXPOTH(u,v)=exp(log(1/sqrt(2))*(D(u,v)/Dcut)^2);

        if D(u,v)<D0

            THPFH(u,v)=1;

        elseif D(u,v)<=D1

            THPEH(u,v)=(D(u,v)-D1)/(D0-D1);

        else

            THPFH(u,v)=0;

        end

    end

end

BUTTERG=BUTTERH.*F;

B=ifft2(BUTTERG);

EXPOTG=EXPOTH.*F;

E=ifft2(EXPOTG);

THPFG=THPFH.*F;

T=ifft2(THPFG);

figure, subplot(221);imshow(noisy); 

title('加噪声图像')%显示图像

subplot(222);imshow(B,map);

 title('巴特沃斯低通滤波')

subplot(223);imshow(E,map) ;

title('指数低通滤波')

subplot(224);imshow(T,map);

 title('梯形低通滤波')