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);%构造滤波函数    endendB=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);%构造滤波函数    endendB=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);%构造滤波函数    endendB=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);    endendE=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);    endendE=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);    endendE=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('梯形低通滤波')