clear;close all;clc;

f=imread('rice.png');
f=im2double(f);

N=[  0 -1  0 ;
    -1  5 -1 ;
     0 -1  0 ]; %锐化模板
[fsha,fedg]=sharpen(f,N); %根据模板进行处理
imshow(fedg);title('已取反的边缘')
figure(2)
subplot(121)
imshow(f);title('原图')
subplot(122)
imshow(fsha);title('锐化图像')

function [fsha,fedg]=sharpen(f,Operator)

%根据模板进行锐化
%适用于真彩色图像与灰度图像
%f为原图像，Operator为锐化模板
%输出的fsha为锐化后的图像，fedg为已取反的加重的边缘
%double型

[m,n,h]=size(f);
fsha=f;

for k=1:h
    for i=2:m-1
        for j=2:n-1
            fsha(i,j,k)=sum(sum(f(i-1:i+1,j-1:j+1,k).*Operator)); %根据模板进行锐化处理
        end
    end
end

if sum(Operator(:))>0
    fedg=f-fsha+1; %取反加重的边缘
else
    fedg=1-fsha; %取反未加重的边缘
    fsha=f+fsha; %此时，锐化图像=原图像+边缘
end

1.2 常用锐化模板

首先我们在MATLAB图库找一张合适的图片，稻米

使用的是

1.2.1 4-邻拉普拉斯锐化模板

N=[0 -1 0 ; -1 5 -1 ; 0 -1 0]

处理后的照片

看一下加重的边缘

进行一下取反，更适合人眼观测

对比

接下来使用

1.2.2 8-邻拉普拉斯锐化模板

看看效果 N=[-1 -1 -1 ; -1 9 -1 ; -1 -1 -1]

很显然，对比度更强了，但不是越强越好，适合人眼观测的才是最好的。

改变锐化强度系数的情况，请自行尝试~

接下来看一下Prewitt算子与Sobel算子(各展示一种)

1.2.3 Prewitt算子

N=[-1 0 1 ; -1 0 1 ; -1 0 1]

1.2.4 Sobel算子

N=[-1 0 1 ; -2 0 2 ; -1 0 1]

相对于Prewitt算子，Sobel算子增强了效果，但对噪声也有放大作用

1.3 水平、垂直、对角方向锐化与平滑

每幅图像都不一样，适宜的锐化模板也不一样，接下来展示单独对某一方向进行锐化处理的结果。

从MATLAB图库选择一张适宜的图片

1.3.1 增强水平方向边缘

N=[0 0 0 ; -1 3 -1 ; 0 0 0]

emmm，加重的边缘可能有些暗淡

那么，加大一下锐化强度系数

N=[0 0 0 ; -2 5 -2 ; 0 0 0]

可以看出，增强水平方向对于水平结构的边缘增强效果一般，但对垂直结构的边缘增强效果良好

1.3.2 增强垂直方向边缘

N=[0 -2 0 ; 0 5 0 ; 0 -2 0]

增强垂直方向则是对水平结构的边缘增强效果较好

1.3.3 增强对角方向边缘

N=[-1 0 -1 ; 0 5 0 ; -1 0 -1]

1.3.4 叠加

锐化强度系数为1时，分别进行水平、垂直、对角方向增强的效果的叠加等效于此8-邻拉普拉斯锐化模板。

1.3.5 增强水平和垂直方向边缘，平滑对角方向边缘

N=[1 -2 1 ; -2 5 -2 ; 1 -2 1]

1.3.6 增强对角方向边缘，平滑水平和垂直方向边缘

N=[-2 1 -2 ; 1 5 1 ; -2 1 -2]

2. 频域(高频提升滤波法)

2.1 几种常用高通滤波器的传递函数

2.2 高通滤波器H(u,v)的特性曲线

截取自《数字图像处理》许录平版 P94 表4.4.1 与 P95 图4.4.6

下面分别进行理想高通滤波(IHPF)，巴特沃思高通滤波(BHPF)及指数高通滤波(EHPF)，三段代码都差不多，改一下传递函数就行

从MATLAB图库里面找一张图

2.3 理想高通滤波(IHPF)

D0可取15，30，50对比一下处理结果，这里仅以D0=15作为展示，下面两种高通滤波也是如此

clear;close all;clc;

f=imread('testpat1.png'); %读取图像
[L,W,H]=size(f);
f=im2double(f); %转换数据类型
F=fftshift(fft2(f)); %将直流分量移到频谱中心
G=zeros(L,W,H);

%理想高通滤波(IHPF)
D0=15; %设置参数
num1=floor(L/2); %向左取整
num2=floor(W/2);
for k=1:H
    for i=1:L
        for j=1:W
            D=sqrt((i-num1)^2+(j-num2)^2);
            if D<=D0
                H=0;
            else
                H=1;
            end
            G(i,j,k)=H*F(i,j,k);
        end
    end
end

F2=ifftshift(G); %去中心化
f2=ifft2(F2); %逆变换
f3=real(f2); %取实数部分
subplot(121);imshow(f);title('原图')
subplot(122);imshow(f3+f);title('理想高通滤波图像');
figure;
subplot(121);imshow(f3);title(sprintf('IHPF滤波结果,D0=%d',D0));
subplot(122);imshow(1-f3);title('取反');

由于理想高通滤波器是突变的，由它得到的高频图像会存在比较强的振铃现象

2.4 巴特沃思高通滤波(BHPF)

clear;close all;clc;

f=imread('testpat1.png'); %读取图像
[L,W,H]=size(f);
f=im2double(f); %转换数据类型
F=fftshift(fft2(f)); %将直流分量移到频谱中心
G=zeros(L,W,H);

%巴特沃斯高通滤波(BHPF)
n=2;D0=15; %设置参数
num1=floor(L/2); %向左取整
num2=floor(W/2);
for k=1:H
    for i=1:L
        for j=1:W
            D=sqrt((i-num1)^2+(j-num2)^2);
            H=1/(1+(2^0.5-1)*(D0/D)^(2*n));
            G(i,j,k)=H*F(i,j,k);
        end
    end
end

F2=ifftshift(G); %去中心化
f2=ifft2(F2); %逆变换
f3=real(f2); %取实数部分
subplot(121);imshow(f);title('原图')
subplot(122);imshow(f3+f);title('巴特沃思滤波图像');
figure;
subplot(121);imshow(f3);title(sprintf('BHPF滤波结果,D0=%d',D0));
subplot(122);imshow(1-f3);title('取反');

巴特沃思高通滤波器的传递函数由于变化比较平滑，在阶数较低时得到的图像只有轻微的振铃现象

2.5 指数高通滤波(EHPF)

clear;close all;clc;

f=imread('testpat1.png'); %读取图像
[L,W,H]=size(f);
f=im2double(f); %转换数据类型
F=fftshift(fft2(f)); %将直流分量移到频谱中心
G=zeros(L,W,H);

%指数高通滤波(EHPF)
n=2;D0=15; %设置参数
num1=floor(L/2); %向左取整
num2=floor(W/2);
for k=1:H
    for i=1:L
        for j=1:W
            D=sqrt((i-num1)^2+(j-num2)^2);
            H=exp(log(1/2^0.5)*(D0/D)^n);
            G(i,j,k)=H*F(i,j,k);
        end
    end
end

F2=ifftshift(G); %去中心化
f2=ifft2(F2); %逆变换
f3=real(f2); %取实数部分
subplot(121);imshow(f);title('原图')
subplot(122);imshow(f3+f);title('指数高通滤波图像');
figure;
subplot(121);imshow(f3);title(sprintf('EHPF滤波结果,D0=%d',D0));
subplot(122);imshow(1-f3);title('取反');