本科毕业论文“基于图像处理的自动对焦技术研究”，对焦过程中的一个重要阶段是图像清晰度评价，我用MATLAB实现了4类清晰度评价函数：基于图像梯度的清晰度评价函数、频域评价函数、信息熵评价函数、统计学评价函数，总计11种。
其中基于图像梯度的清晰度评价函数我一部分借鉴他人上传的毕设代码，在此附上链接。
http://www.pudn.com/Download/item/id/1394579.html
一、基于图像梯度的清晰度评价函数
这类算法的实质是图像处理知识
正焦的清晰图像比模糊的离焦图像边缘更加锐利清晰，边缘像素灰度值变化大，因而有更大的梯度值。在进行图像处理时，将图像看作二维离散矩阵，利用梯度函数获取图像灰度信息，以此来评判图像清晰度。在离散信号中梯度表现为差分形式。常用的梯度函数有：能量梯度函数EOG、Roberts函数、Tenengrad函数、Brenner函数、方差Variance函数、拉普拉斯Laplace函数等等。下面分别对其进行介绍，其中表示对应像素点的灰度值。
1.能量梯度函数（Energy of Gradient,EOG）
将 x 方向和 y 方向的相邻像素的灰度值之差[]的平方和作为每个像素点的梯度值，对所有像素梯度值累加作为清晰度评价函数值，表达式如下所示：
matlab代码如下：
%EOG(Energy Of Grad)N1 = 5; %要处理的图片张数A = zeros(1,N1); %存储每一幅图像清晰度评价值X = zeros(1,N1); %存储做归一化处理后的评价值
tic %计时
for L=1: N1 I=imread([int2str(L),'.jpg']); %连续读取图片I=double(I); [M N]=size(I); FI=0; for x=1:M-1 for y=1:N-1 % x方向和y方向的相邻像素灰度值只差的的平方和作为清晰度值FI=FI+(I(x+1,y)-I(x,y))*(I(x+1,y)-I(x,y))+(I(x,y+1)-I(x,y))*(I(x,y+1)-I(x,y));end end A(1,L) = FI;
end
time=toc
%--------------------------------
%对图像清晰度值做归一化处理，线性函数归一化公式for W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end
%做曲线拟合输出函数曲线
x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2);
Y=polyconf(p,x1,y1);
plot(x1,y1,'y');
hold on;2.Roberts函数
Roberts函数与能量梯度函数相似，它是利用对角方向像素点灰度值之差。将4个相邻像素点的灰度值交叉相减的平方和作为每个像素点的梯度值，对所有像素梯度值累加作为清晰度评价函数值，表达式如下式所示：
%RobertsN1 = 5; A = zeros(1,N1); X = zeros(1,N1);ticfor L=1: N1 I=imread([int2str(L),'.jpg']); I=double(I); [M N]=size(I); FI=0; %Robert算子原理，对角方向相邻的两像素之差 for x=1:M-1 for y=1:N-1 FI= FI + (abs(I(x,y)-I(x+1,y+1))+abs(I(x+1,y)-I(x,y+1))); end end A(1,L) = FI; end time=tocfor W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2);
Y=polyconf(p,x1,y1);
plot(x1,y1,'c');
hold on;3.Tenengrad函数
采用Sobel算子提取像素点水平方向和垂直方向的梯度值，Tenengrad函数定义为像素点梯度的平方和，并为梯度设置一个阈值T调节函数的灵敏度。表达式如下式所示：
代码：
%TenengradN1 =5; A = zeros(1,N1); X = zeros(1,N1); ticfor L=1: N1 I=imread([int2str(L),'.jpg']); I=double(I); [M N]=size(I); %利用sobel算子gx,gy与图像做卷积，提取图像水平方向和垂直方向的梯度值
GX = 0; %图像水平方向梯度值
GY = 0; %图像垂直方向梯度值
FI = 0; %变量，暂时存储图像清晰度值
T = 0; %设置的阈值for x=2:M-1 for y=2:N-1 GX = I(x-1,y+1)+2*I(x,y+1)+I(x+1,y+1)-I(x-1,y-1)-2*I(x,y-1)-I(x+1,y-1); GY = I(x+1,y-1)+2*I(x+1,y)+I(x+1,y+1)-I(x-1,y-1)-2*I(x-1,y)-I(x-1,y+1); SXY= sqrt(GX*GX+GY*GY); %某一点的梯度值%某一像素点梯度值大于设定的阈值，将该像素点考虑，消除噪声影响if SXY>T FI = FI + SXY*SXY; %Tenengrad值定义end end end A(1,L) = FI; end time=toc% X = zeros(1,N1); for W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end
x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2);
Y=polyconf(p,x1,y1);
plot(x1,y1,'g');
hold on;
4.Brenner函数
又称梯度滤波器法，只需计算x方向上相差两个像素点的差分，即计算二阶梯度，计算量少，表达式如下所示：
代码：
%BrennerN1 = 5; %N1为要处理的图片张数A = zeros(1,N1); %zeros()定义指定行列的零矩阵；A矩阵用来存储每一幅图像的清晰度原值X = zeros(1,N1); %X用来存储做归一化处理后的函数值%------------------------------ticfor L=1: N1 I=imread([int2str(L),'.jpg']); %读取图片，将值转换为字符串接受向量和矩阵输入
I=double(I); %精度存储问题[M N]=size(I); %M等于矩阵行数，N等于矩阵列数；size()获取矩阵行列FI=0; %变量，暂时存储每一幅图像的Brenner值for x=1:M-2 %Brenner函数原理，计算相差两个位置的像素点的灰度值for y=1:N FI=FI+(I(x+2,y)-I(x,y))*(I(x+2,y)-I(x,y)); end end A(1,L) = FI; endtime=toc%对原始数据做归一化处理，线性函数归一化公式for W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end
%曲线拟合
x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2); %polyfit(x,y,n)曲线拟合函数，已知离散点坐标拟合曲线；x,y横纵坐标，n为拟合阶数，一阶直线拟合，二阶抛物线拟合 ，返回幂次从高到低的多项式系数向量P，矩阵S用于生成预测值的误差估计
Y=polyconf(p,x1,y1); %置信区间
plot(x1,y1,'r'); %画出拟合曲线，红线red
title('梯度评价函数');
xlabel('成像面位置');
ylabel('归一化后的图像清晰度评价值');
hold on;5.Variance方差函数
Variance函数表示图像灰度分布的离散程度。离焦图像灰度值变换范围小，离散程度低，方差小；正焦图像灰度值变换范围大，离散程度高，方差大。因此可以用其作为评价函数，对于M*N大小的图像表达式如下式所示：
代码：
%VarianceN1 = 5; A = zeros(1,N1); X = zeros(1,N1);
ticfor L=1: N1I=imread([int2str(L),'.jpg']); I=double(I); [M N]=size(I); gama = 0; %gama图像平均灰度值%求gamafor x=1:M for y=1:N gama = gama + I(x,y); end end gama = gama/(M*N); FI=0; for x=1:M for y=1:N FI=FI+(I(x,y)-gama)*(I(x,y)-gama); end end A(1,L) = FI;end time=tocfor W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];[p,S]=polyfit(x1,y1,2); Y=polyconf(p,x1,y1); plot(x1,y1,'b');hold on;6. Laplace拉普拉斯函数
采用Laplace算子与图像各个像素点的灰度值进行卷积得到一个梯度矩阵记为，取各像素点梯度的平方和作为评价函数，如下式所示：
代码：
%Laplace N1 = 5; A = zeros(1,N1); X = zeros(1,N1);ticfor L=1: N1 I=imread([int2str(L),'.jpg']); I=double(I); [M N]=size(I); FI=0; for x=2:M-1 for y=2:N-1 IXXIYY = -4*I(x,y)+I(x,y+1)+I(x,y-1)+I(x+1,y)+I(x-1,y); FI=FI+IXXIYY*IXXIYY; %取各像素点梯度的平方和作为清晰度值 end end
A(1,L) = FI; end time=tocfor W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2);
Y=polyconf(p,x1,y1);
plot(x1,y1,'m');
hold off;二：频域评价函数
1.基于二维离散傅里叶变换的图像清晰度评价函数
其中，加权系数表示像素到中心像素的距离，其作用为强调频谱中的高频成分。因为图像经过傅里叶变换和移位后，频率分布为从中央到四周边缘从低频向高频扩散；M,N 为图像尺寸；P（u，v）表示图像的频谱的平方。
代码：
%DFT
N1 = 5; %N1要处理的图像张数
A = zeros(1,N1); %A向量用来存储每一幅图像的清晰度原值
X = zeros(1,N1); %X用来存储做归一化处理后的清晰度函数值
%----------------------
tic
for L=1: N1
I=imread([int2str(L),'.jpg']);
I=rgb2gray(I);
I=double(I);
[M N]=size(I);
fftI = fft2(I); %进行二维离散傅里叶变换
sfftI = fftshift(fftI); %移位，直流分量移到图像中心
magnitude = abs(sfftI); %取模值
FI=0;
for u=1:Mfor v=1:NFI=FI+sqrt(u*u+v*v)*magnitude(u,v); %基于离散傅里叶变换的清晰度评价函数end
end
A(1,L) = FI/(M*N);
end
time=toc
%对原始数据做归一化处理，线性函数归一化公式
%-------------------------for W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end
%曲线拟合，0为正焦位置，离焦-正焦-离焦
x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2); %polyfit(x,y,n)曲线拟合函数，已知离散点坐标拟合曲线；x,y横纵坐标，n为拟合阶数，一阶直线拟合，二阶抛物线拟合 ，返回幂次从高到低的多项式系数向量P，矩阵S用于生成预测值的误差估计
Y=polyconf(p,x1,y1); %置信区间
plot(x1,y1,'b'); %画出拟合曲线，红线red
title('频域评价函数');
xlabel('成像面位置');
ylabel('归一化后的图像清晰度评价值');
hold on;2.离散余弦变换DCT
基于DFT的评价函数灵敏度高，但傅里叶变换系数是复数，计算量大。离散余弦变换DCT对此进行改进，变换系数是实数，计算量减少且同样可以表示图像频率分布信息。对于一幅M*N的图像，基于离散余弦变换DCT变换的清晰度评价函数定义如下：
%DCT
N1 = 5; %N1要处理的图像张数
A = zeros(1,N1); %A向量用来存储每一幅图像的清晰度原值
X = zeros(1,N1); %X用来存储做归一化处理后的清晰度函数值
%----------------------
tic
for L=1: N1
I=imread([int2str(L),'.jpg']);
I=rgb2gray(I);
I=double(I)+10*randn(size(I));
[M N]=size(I);
dctI = dct2(I); %进行二维离散傅里叶变换
magnitude = abs(dctI); %取模值
FI=0;
for u=1:Mfor v=1:NFI=FI+(u+v)*magnitude(u,v); %基于离散傅里叶变换的清晰度评价函数end
end
A(1,L) = FI/(M*N);
end
time=toc
%对原始数据做归一化处理，线性函数归一化公式
%-------------------------for W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end
%曲线拟合，0为正焦位置，离焦-正焦-离焦
x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2); %polyfit(x,y,n)曲线拟合函数，已知离散点坐标拟合曲线；x,y横纵坐标，n为拟合阶数，一阶直线拟合，二阶抛物线拟合 ，返回幂次从高到低的多项式系数向量P，矩阵S用于生成预测值的误差估计
Y=polyconf(p,x1,y1); %置信区间
plot(x1,y1,'r'); %画出拟合曲线，红线red
hold off;
三、基于信息熵的清晰度评价函数
在信息论中，熵用来描述信息的丰富程度。信息熵评价函数基于正焦图像中具有多样性的灰度分布，像素灰度值分布区间广，灰度值之间差异大，此时熵值大；离焦图像正好相反。因此可以基于信息熵构造清晰度评价函数，定义如下：
代码：
%entropy
N1 = 5; %N1为要处理的图片张数
A = zeros(1,N1); %zeros()定义指定行列的零矩阵；A矩阵用来存储每一幅图像的清晰度原值
X = zeros(1,N1); %X用来存储做归一化处理后的函数值
%处理图片
tic
for L=1: N1 I=imread([int2str(L),'.jpg']); %读取图片，将值转换为字符串接受向量和矩阵输入I=rgb2gray(I);I=double(I); A(1,L)=entr(I); %调用求熵值函数
end
time=toc
%对原始数据做归一化处理，线性函数归一化公式for W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end
x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2); %polyfit(x,y,n)曲线拟合函数，已知离散点坐标拟合曲线；x,y横纵坐标，n为拟合阶数，一阶直线拟合，二阶抛物线拟合 ，返回幂次从高到低的多项式系数向量P，矩阵S用于生成预测值的误差估计
Y=polyconf(p,x1,y1); %置信区间
plot(x1,y1,'r'); %画出拟合曲线，红线red
title('基于信息熵的评价函数');
xlabel('成像面位置');
ylabel('归一化后的图像清晰度评价值');%定义子函数entr()，求一幅图像的熵值
%-------------------------------------------
function[H_img]= entr(I)
[C,R]=size(I); %求图像的规格
Img_size=C*R; %图像像素点的总个数
L=256; %图像的灰度级0-255
H_img=0; %图象熵
nk=zeros(L,1); %存储图像灰度出现次数
for i=1:C
for j=1:R
Img_level=I(i,j)+1; %获取图像的灰度级
nk(Img_level)=nk(Img_level)+1; %统计每个灰度级像素的点数
end
end
for k=1:L
Ps(k)=nk(k)/Img_size; %计算每一个灰度级像素点所占的概率
if Ps(k)~=0 %去掉概率为0的像素点
H_img=-Ps(k)*log2(Ps(k))+H_img; %求熵值的公式
end
end
end四、基于统计学的清晰度评价函数
1.Range函数
Range函数基于灰度直方图中灰度带的变化大小反映清晰度，定义如下所示：
代码：
%Range
N1 =5;
gray_level = 32; %灰度直方图中划分的灰度等级
temp=zeros(1,gray_level);
A = zeros(1,N1);
X = zeros(1,N1);
tic
for L=1: N1 I=imread([int2str(L),'.jpg']); %读取图片，将值转换为字符串接受向量和矩阵输入I=rgb2gray(I); I=double(I); %-------没做数值类型，出错，曲线相反----------[count,K] = imhist(I,gray_level);%imhist()画灰度分布直方图,count表示某一灰度区间的像素个数，K表示灰度区间取值for y=1:gray_level temp(1,y)=count(y)*K(y);endA(1,L)=max(temp)-min(temp);
end
time=toc
for W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R;
end
x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2); %polyfit(x,y,n)曲线拟合函数，已知离散点坐标拟合曲线；x,y横纵坐标，n为拟合阶数，一阶直线拟合，二阶抛物线拟合 ，返回幂次从高到低的多项式系数向量P，矩阵S用于生成预测值的误差估计
Y=polyconf(p,x1,y1); %置信区间
plot(x1,y1,'g'); %画出拟合曲线，红线red
title('统计学评价函数');
xlabel('成像面位置');
ylabel('归一化后的图像清晰度评价值');
hold on2 Vollaths函数
Vollaths函数又称自相关函数，反映空间两点的相似性。正焦图像边缘清晰锐利，像素点之间相关程度低；离焦图像像素点相关程度高。清晰度评价函数如下所示：
.代码：
%vollathsN1 = 5; %N1为要处理的图片张数A = zeros(1,N1); %zeros()定义指定行列的零矩阵；A矩阵用来存储每一幅图像的清晰度原值X = zeros(1,N1); %X用来存储做归一化处理后的函数值%用一个for循环处理每一张图片ticfor L=1: N1 I=imread([int2str(L),'.jpg']); %读取图片，将值转换为字符串接受向量和矩阵输入I=double(I); %精度存储问题[M N]=size(I); %M等于矩阵行数，N等于矩阵列数；size()获取矩阵行列数%begintime=clock; FI=0; %变量，暂时存储每一幅图像的Brenner值for x=1:M-2 %Brenner函数原理，计算相差两个位置的像素点的灰度值for y=1:N FI=FI+I(x,y)*abs(I(x+1,y)-I(x+2,y)); end end %time=etime(clock,begintime); A(1,L) = FI; end time=toc%对原始数据做归一化处理，线性函数归一化公式for W = 1:N1 C = max(A); D = min(A); E = C-D; R = (A(1,W) - D)/(E); X(1,W) = R; end x1=[-20 -10 0 10 20 ];
y1 = [X(1,1) X(1,2) X(1,3) X(1,4) X(1,5)];
[p,S]=polyfit(x1,y1,2); %polyfit(x,y,n)曲线拟合函数，已知离散点坐标拟合曲线；x,y横纵坐标，n为拟合阶数，一阶直线拟合，二阶抛物线拟合 ，返回幂次从高到低的多项式系数向量P，矩阵S用于生成预测值的误差估计
Y=polyconf(p,x1,y1); %置信区间
plot(x1,y1,'r'); %画出拟合曲线，红线red
hold off以下是博主的仿真结果：
梯度
频域
信息熵
统计学
本科毕业论文，完全小白，不足之处，敬请谅解。
