(完整版)DCT变换与KLT变换在图像压缩中的应用

DCT变换与KLT变换在图像压缩中的应用
专业：电子与通信工程
学号：16213026 16213024
姓名：
目录
1 报告简介 (1)
2 算法原理 (1)
2.1 DCT算法原理 (1)
2.2 KLT算法原理 (3)
3 仿真分析 (5)
3.1 DCT仿真分析 (5)
3.2 KLT仿真分析 (7)
5 总结 (8)
参考文献 (9)
附录：MATLAB代码 (10)
DCT变换与KLT变换在图像压缩中的应用
1 报告简介
随着数据化时代的开启，图像压缩技术越来越成熟并且应用越来越广泛，本文在研究近年来图像压缩一般方法的基础上，介绍了基于DCT变换的图像压缩的基本原理及其实现步骤，以及扩展研究了KL变换的图像压缩方法，并使用MATLAB,针对同一幅原始图像进行不同方法的压缩比较，给出了实验仿真结果。

本论文首先提出了用MATLAB来实现DCT变换的数字图像压缩技术, 方法简单,快速,且误差小。

然后介绍了KLT图像压缩编码的具体过程和方法。

最后分析了图像经过2种压缩方法时，图像质量的变化情况。

2 算法原理
2.1 DCT算法原理
DCT变换利用傅立叶变换的性质。

采用图像边界褶翻将像变换为偶函数形式，然后对图像进行二维傅立叶变换，变换后仅包含余弦项，所以称之为离散余弦变换。

DCT编码属于正交变换编码方式，用于去除图像数据的空间冗余。

变换编码就是将图像光强矩阵(时域信号)变换到系数空间(频域信号)上进行处理的方法。

在空间上具有强相关的信号，反映在频域上是在某些特定的区域内能量常常被集中在一起，或者是系数矩阵的分布具有某些规律。

我们可以利用这些规律在频域上减少量化比特数，达到压缩的目的。

图像经DCT变换以后，DCT系数之间的相关性就会变小。

而且大部分能量集中在少数的系数上，因此，DCT变换在图像压缩中非常有用，是有损图像压缩国际标准JPEG的核心。

从原理上讲可以对整幅图像进行DCT变换，但由于图像各部位上细节的丰富程度不同，这种整体处理的方式效果不好。

为此，发送者首先将输入图像分解为8
8⨯或16
16⨯块，然后再对每个图像块进行二维DCT 变换，接着再对DCT系数进行量化、编码和传输；接收者通过对量化的DCT系数
进行解码，并对每个图像块进行的二维DCT 反变换。

最后将操作完成后所有的块拼接起来构成一幅单一的图像。

对于一般的图像而言，大多数DCT 系数值都接近于0，所以去掉这些系数不会对重建图像的质量产生较大影响。

因此，利用DCT 进行图像压缩确实可以节约大量的存储空间。

基于DCT 的JPEG 图像压缩编码算法原理可用图1表示：
（a ）编码器
（b ）解码器 图1 DCT 算法原理框图
在编码过程中，首先将输入图像颜色空间转换后分解为8×8大小的数据块，然后用正向二维
DCT 把每个块转变成64个DCT 系数值，其中1个数值是直流(DC)系数，即8×8空域图像子块的平均值，其余的63个是交流(AC)系数，接下来对DCT 系数进行量化，量化过程实际上就是对 DCT 系数的一个优化过程，它是利用了人眼对高频部分不敏感的特性来实现数据的大幅简化。

量化过程实际上是简单地把频率领域上每个成份，除以一个对于该成份的常数，且接着四舍五入取最接近的整数。

这是整个过程中的主要有损运算。

量化是图像质量下降的最主要原因。

量化后的数据，有一个很大的特点，就是直流分量相对于交流分量来说要大，而且交流分量中含有大量的0。

这样，对这个量化后的数据如何来进行简化，从而再更大程度地进行压缩呢。

将量化后的系数进行“Z ”字形编排，这样做的特点就是会连续出现多个0，即充分利用相邻两图像块的特性，来再次简化数据，从而再更大程度地进行压缩。

最后将变换得到的量化的DCT 系数进行编码和传送，这样就完成了图像的压缩过程。

压缩图像
在解码过程中，形成压缩后的图像格式，先对已编码的量子化的DCT 系数进行解码，然后求逆量化并把DCT 系数转化为8×8样本像块(使用二维DCT 反变换)，最后将操作完成后的块组合成一个单一的图像。

这样就完成了图像的解压过程。

一个N N ⨯块),(y x f 的二维DCT ),(v u F 的定义如下：
∑∑-=-=++=1010
)2)12(cos(),()2)12(cos()()(2),(N x N y N v y y x f N u x v C u C N v u F ππ （1） 对应的N N ⨯块的二维IDCT 则为：
∑∑-=-=++=1010)2)12(cos(),()2)12(cos(
)()(2),(N u N v N v
y v u F N u x v C u C N y x f ππ （2）
式中，空域的x 、y ，频域的u 、v 取值集合均为}1,1,0{-N ，其中：
⎪⎩⎪⎨⎧==其他10,u 21)(),(v v C u C （3）
将离散余弦变换变换写为矩阵形式为：
T N N fC C F = （4）
N T N FC C f = （5）
其中，N C 为N N ⨯正交变换矩阵，f 为N N ⨯原图像块，F 为N N ⨯变换域图像块。

2.2 KLT 算法原理
KLT 变换又称Hotelling 变换，特征向量变换或主分量方法。

KL 变换是遥感图像增强和信息提取中用得最多的线性变换，是对原波段图像进行波谱信息的线性投影变换，在尽可能不减少信息量的前提下，将原图像的高维多光谱空间的像元亮度值投影到新的低维空间，减少特征空间维数，达到数据压缩、提高信噪比、提取相关信息、降维处理和提取原图像特征信息的目的，并能有效地提取影像信息。

它可使原来多波段图像经变换后提供出一组不相关的图像变量，最前面的主分量具有较大的方差，包含了原始影像的主要信息，所以要集中表达信息，突出图像的某些细部特征，可采用主分量变换来完成。

KL 变换是图像压缩中的一种最优正交变换。

KLT 的突出优点是去相关性好，它根据具体的图像统计特性来决定它的变换矩阵，对图像有最好的匹配效果，能将信号在变换域的相关性全部解除，是最小均方误差下的最佳变换。

图2 KLT 算法原理框图
KLT 变换就是对8x8的图像矩阵求自协方差矩阵，对自协方差矩阵做特征值分解，得到由特征值从小到大排列的对角矩阵，和由特征向量组成的矩阵，特征矩阵与图像矩阵做乘法，称为正交变换，即KL 变换，的到的新的矩阵每一行称为一个新的变量，其中第一行几乎包含了总方差%80以上的信息，其余行包含的信息依次减少，新矩阵的个元素之间是不相关的，因而KL 变换去掉了变量之间的相关性。

KLT 是对向量x 做的一个正交变换x y Φ=，目的是变换到y 后去除数据相关性。

其中，Φ是x 特征向量组成的矩阵，满足I T =ΦΦ，当x 都是实数时，Φ是正交矩阵。

用y m 表示向量y 的平均值，y 的协方差矩阵记为∑y ，通过变换Tx y Φ=，得到：
),,,(])([])([))((]))([()(110-=Λ=ΦΦ=Φ-Φ=ΦΦ-ΦΦ=ΦΦ-ΦΦ=-=∑∑
N T T
x x T T T x x T T T T x T x T T T T T y diag E E E E λλλ x T y y m m xx m m xx m m x x m m yy （6）
写成矩阵形式：
⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=Λ=ΦΦ=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=--∑∑2121201100000000000000N N T ij y σσσλλλσ x （6）
由此可见，做了KLT 变换之后呢，∑y 成为了对角阵，也就是对于任意j i ≠，有0),cov(=j i y y ，当j i =有i j i y y λ=),cov(，因此利用KLT 去除了数据相关性。

而且，i y 的方差与x 协方差矩阵的第i 个特征值相等即i i λσ=2。

3 仿真分析
利用MATLAB 对算法进行仿真分析，测试环境为：操作系统为Win7，CPU 为i5-3210M ，内存为4GB ，使用的MATLAB 版本为R2012b 。

利用MATLAB 产生仿真图形界面如图3所示。

图3 DCT 和KLT 仿真图形界面设计
3.1 DCT 仿真分析
输入一幅图像，通过改变DCT 量化值进行图像压缩，同时，利用相关参数对图
像质量进行对比分析，图4为经过DCT压缩的效果对比图。

（a）原始图像（b）量化程度为16
（c）量化程度为32 （b）量化程度为64
图3 运用DCT变化进行图像压缩效果对比图
通过肉眼观察，可以发现量化程度越大图像越清晰，反之，越模糊，细节信息越差。

采用不同的量化程度可以改变图像的清晰度及文件的大小，本文将利用图像的相关性能参数对压缩图像进行比较。

图像的方差、平均梯度、信息熵及对比度与图像质量正相关。

图像方差越大，表明图像灰度层次越丰富；平均梯度越大，图像层次越多；信息熵越大，表明图像的信息量越大。

通过表1，可以清晰地看出量化程度越大，得到的压缩图像性能越好，相对的图像大小也增大，但原图相比，图像性能变化不明显，但压缩效果明显。

表1 DCT变换得到压缩图像性能参数比较
方差平均梯度信息熵（bit）文件大小（KB）图3（a） 1.9545e+03 5.4304 7.1771 26.4
图3（b） 1.8785e+03 3.4470 7.1412 9.82
图3（c） 1.9011e+03 4.1021 7.1538 12
图3（d） 1.9167e+03 4.3852 7.1558 12.7
3.2 KLT仿真分析
利用KLT算法对图4（a）图进行图像压缩得到图5。

（a）量化程度为16 （b）量化程度为32
（c）量化程度为32 （b）量化程度为64
图4运用KLT变化进行图像压缩效果对比图
从图像的方差、信息熵及平均梯度对图像性能进行分析，得到表2。

表2 KLT变换得到压缩图像性能参数比较
由表2可以得到KLT变换可以对图像进行压缩，并且采用相同的量化程度，DCT 比KLT变化效果较好。

5 总结
本文叙述的图像压缩编码算法，DCT是正交变换，它可以将8*8图像的空间表达式转换为频率域，只需要用少量的数据点表示图像；DCT产生的系数很容易被量化，因此能获得好的块压缩；DCT算法的性能很好，它有快速算法，因此它在软件中容易实现；而且DCT算法是对称的，所以利用逆DCT算法可以用来解压缩图像。

由实验仿真结果可知，当图像压缩比增大时，也即压缩效率减小时，图像的质量也将降低，人们可以根据需要的图像的质量来规定压缩比的大小。

图像经过不同的压缩比后，图像的质量变化的同时，图像的信噪比也跟着变化。

压缩比增大时，则在信道传输的时候丢失的信息就越多，这样使得信号与噪声的比值变小。

利用DCT 变换进行图像压缩可以节约大量的存储空间。

压缩应该在最合理地近似原图像的情况下使用最少的系数。

由于之前存在对KLT变换的误解，在进理论补充后，发现该方法并不能全部消除变换后系数间的相关性，各数据块之间相关性没有消除，数据块也可能存在极弱的相关性。

由于水平有限，并不能解决，即使如此，该方法应用于图像的压缩编码中效果也极其明显。

KLT没有快速算法，因此对信息量大的图像编码速度比较慢，但在此次试验中用的是较小的图像，所以处理速度效果不明显。

参考文献
[1]N Ahmed，T Natarajan，K.T.Rao．Discrete cosine transform[J]．IEEE Tran．On Computers, 1974, C-23: 90-93．
[2] 全子一．图像信源压缩编码及信道传输理论与新技术[M]．北京：北京工业大学出版社，
2006．
[3] 何小海，主编．数字图像通信及其应用[M]．成都:四川大学出版社，2006．
[4] 罗开仲，黄士坦，杨华民．DCT算法及其与小波编码在图像处理中的比较[J]．计算机技术
与发展，2006，16(9)．
[5] Xiong Z，Orchard M，Guleryuz O．A DCT-based embedded image coder [J]．IEEE Signal Processing Letters，1996，11(3):289-290．
[6]Shaorong Chang，Lawrence Carin，A modified SPIHT algorithm for image coding with
a joint MSE and classification distortion measure [J]. IEEE Transactions on Image Processing，2006，15(3): 713-725．
附录：MATLAB代码
%%主程序
clc;
clear
I=imread('football.jpg');
clf reset
set(gcf,'unit','normalized','position',[0.1,0.2,0.45,0.45]);%设置图形窗大小
set(gcf,'defaultuicontrolunits','normalized');
set(gcf,'defaultuicontrolfontsize',11);
set(gcf,'defaultuicontrolfontname','楷书');
set(gcf,'defaultuicontrolhorizontal','left');
% set(gcf,'menubar','none'); %删除图形窗工具条str='KLT DCT性能比较';
set(gcf,'name',str ,'numbertitle','off'); %书写图形窗名
z=32;
str1='量化程度';
htext=uicontrol(gcf,'style','text',...%制作静态说明文本框
'position',[0.1,0.8,0.33,0.1],...
'string',[str1,sprintf('%1.4g\',z)]);
hslider=uicontrol(gcf,'style','slider',...%创建滑动键
'position',[0.1,0.6,0.8,0.1],...
'max',64,'min',0,...%设最大阻尼比为2，最小阻尼比为0.02
'sliderstep',[1/64,1/64],...%箭头操纵滑动步长1%，游标滑动步长5%
'Value',32); %缺省取阻尼比等于0.5
set(hslider,'callback',[...%操作滑动键,引起回调'z=get(gcbo,''value'');',...%获得滑动键状态值
'callcheck(htext,str1,z)']); %被回调的函数文件
hpush=uicontrol(gcf,'Style','push',...%制作与列表框配用的按键
'position',[0.1,0.05,0.18,0.15],'string','原始图像');
set(hpush,'callback','jpeg_show(I)')
hpush=uicontrol(gcf,'Style','push',...%制作与列表框配用的按键
'position',[0.4,0.05,0.18,0.15],'string','KLT');
set(hpush,'callback','jpeg_klt(I,z)')
hpush=uicontrol(gcf,'Style','push',...%制作与列表框配用的按键'position',[0.7,0.05,0.18,0.15],'string','DCT');
set(hpush,'callback','jpeg_dct(I,z)')
%灰度图像的klt正反变换
function y=klt(x,mask)
IM=x;%读取图像
IM=im2uint8(IM);
%若图像尺寸不是8的整数倍，补0，补成8的整数倍
[cc1,cc2]=size(IM);%cc1为IM的行数，cc2为IM的列数
a=8-mod(cc1,8);
b=8-mod(cc2,8);
cc3=cc1+a;
cc4=cc2+b;
IM(cc3,cc4)=IM(cc1,cc2);
IM(1:cc1,cc2+1:cc4)=IM(1:cc1,cc2-b+1:cc2);
IM(cc1+1:cc3,1:cc2)=IM(cc1-a+1:cc1,1:cc2);
IM(cc1+1:cc3,cc2+1:cc4)=IM(cc1,cc2);
%若图像尺寸不是8的整数倍，补成8的整数倍
k=1;
for l1=1:8:(cc3-7) %将图像矩阵--8X8的矩阵存起来for l2=1:8:(cc4-7)
F(:,:,k)=IM(l1:l1+7,l2:l2+7);
k=k+1;
end%将图像矩阵--8X8的矩阵存起来end
F=uint8(F);
k=k-1;
qq=k;
d=0;
d=uint32(d);
IM=uint32(IM);
for n1=1:cc3 %求均值
for n2=1:cc4
d=IM(n1,n2)+d;
end
end
d=d/(cc1*cc2);
d=uint8(d); %求均值
IM=uint8(IM);
B2=zeros(64);
jz=d;%jz为灰度均值
for ee=1:qq
A=F(:,:,ee);
TZZ;%调用函数求一个8X8的矩阵的特征矩阵
B2=B2+B;
end
B=B2/ee;%B为均方差矩阵
[V,D]=eig(B);%求特征向量与特征值
Q=V';
for k=1:qq
A=F(:,:,k);
ee=1;
for n1=1:8 %将元素变换成一维
for n2=1:8 %将元素变换成一维
X(ee,1)=A(n1,n2);
YSZ(ee,1)=mask(n1,n2);%YSZ为mask变换成的一维数组
ee=ee+1;
end%将元素变换成一维
end%将元素变换成一维
X=double(X);
Y=Q*X;%变换后的矩阵放在Y中
Y=flipud(Y);
Y=Y.*YSZ;
Y=flipud(Y);
X1(:,:,k)=Q'*Y;%求还原后的矩阵
end
k=1;
for l1=1:8:(cc3-7) %将矩阵8X8矩阵放入IM2中
for l2=1:8:(cc4-7)%将矩阵8X8矩阵放入IM2中
H2=X1(:,:,k);
ee=1;
for n1=1:8
for n2=1:8
C2(n1,n2)=H2(ee,1);
ee=ee+1;
end
end
IM2(l1:l1+7,l2:l2+7)=C2;
k=k+1;
end%将矩阵8X8矩阵放入IM2中
end%将矩阵8X8矩阵放入IM2中
IM=IM(1:cc1,1:cc2);
IM2=IM2(1:cc1,1:cc2);
y=uint8(IM2);
function jpeg_dct(I,z)
tic
% I=imread('football.jpg');%读入原始图像，该图片在安装matlab的目录中找，原图为jpeg图象
Y=im2double(I);%图像存储类型转换，matlab读入图像的数据是uint8，而matlab中数值一般采用double 型（64位）存储和运算。

所以要先将图像转为double格式的才能运算
T=dctmtx(8);%离散余弦变换矩阵，使用由dctmtx函数返回的DCT变换矩阵，这种方法较适合于较小的输入方阵（例如8×8或16×16），计算二维DCT，矩阵T及其转置是DCT函数P1*X*P2的参数%B = blkproc(A,[m n],fun, parameter1, parameter2, ...)
%x就是每一个分成的8*8大小的块，P1*x*P2相当于像素块的处理函数，p1=T p2=T’,也就是
fun=p1*x*p2'=T*x*T'的功能是进行离散余弦变换
%x代表A，对一副原始图像，选取它的最左上角8x8的矩阵，并将每一个像素值转换到-128到127的范围内，得到矩阵x.
%RGB图像分三个通道，分别处理
B(:,:,1)=blkproc(Y(:,:,1),[8 8],'P1*x*P2',T,T');%将原始图像8x8的像素块转换成代表不同频率分量的系数集，DCT后的64个DCT系数与DCT前的64个像素块相对应
B(:,:,2)=blkproc(Y(:,:,2),[8 8],'P1*x*P2',T,T');
B(:,:,3)=blkproc(Y(:,:,3),[8 8],'P1*x*P2',T,T');
%对原图像进行DCT变换
z=fix(z);%量化，向0靠拢取整
t(1:z)=1;
t(z+1:64)=0;
mask=jpeg_Dzz8(t);%Z字形扫描
B2(:,:,1)=blkproc(B(:,:,1),[8 8],'P1.*x',mask);%只保留DCT变换的左上角10个系数，数据压缩，丢弃右下角高频数据，达到图像压缩的目的
B2(:,:,2)=blkproc(B(:,:,2),[8 8],'P1.*x',mask);%变换后的系数值较大的会集中在区域的左上部，即低频分量都集中在左上部。

保留的也是这一部分。

B2(:,:,3)=blkproc(B(:,:,3),[8 8],'P1.*x',mask);%其他部分的系数被舍去，在恢复信号时对它们补0 。

这样以来，由于保留了大部分图像信号能量，在恢复信号后，其质量不会产生显著变化。

%数据压缩，丢弃右下角高频数据
I2(:,:,1)=blkproc(B2(:,:,1),[8 8],'P1*x*P2',T',T);%重构图像
I2(:,:,2)=blkproc(B2(:,:,2),[8 8],'P1*x*P2',T',T);
I2(:,:,3)=blkproc(B2(:,:,3),[8 8],'P1*x*P2',T',T);
%进行DCT反变换，得到压缩后的图像
figure;
set(gcf,'unit','normalized','position',[0.4,0.5,0.3,0.3]);%设置图形窗大小
set(gcf,'menubar','none'); %删除图形窗工具条
str='量化程度';
set(gcf,'name',[str sprintf('%1.4g\',z)],'numbertitle','off'); %书写图形窗名
imshow(I2)
title('DCT变换压缩后的图像')
imwrite(I2,'dct.jpg');%输出压缩后的图像，文件名为football_dct.jpg
a=toc;
uicontrol('style', 'text', 'unit', 'normalized', ...
'position', [0 0 0.2 0.1], 'string', num2str(a), ...
'FontSize', 18)
%搜索方式：对角的Z字形扫描;二维的图像矩阵利用一维的搜索方式。

function y=jpeg_Dzz8(x)
t=1;
for k=2:9
for m=1:k-1
if rem(k,2)==0
y(k-m,m)=x(t);
t=t+1;
else
y(m,k-m)=x(t);
t=t+1;
end
end
end
for k=10:16
for m=k-8:8
if rem(k,2)==0
y(k-m,m)=x(t);
t=t+1;
else
y(m,k-m)=x(t);
t=t+1;
end
end
end
function callcheck(htext,str1,z)
set(htext,'string',[str1,sprintf('%1.4g\',z)]); %更新静态文本框内容<2
function jpeg_klt(I,z)
tic
% I=imread('football.jpg');%读入原始图像，该图片在安装matlab的目录中找，原图为jpeg图象z=fix(z);%量化，，向0靠拢取整
t(1:z)=1;
t(z+1:64)=0;
mask=jpeg_Dzz8(t);%Z字形扫描
I2(:,:,1)=klt(I(:,:,1),mask);
I2(:,:,2)=klt(I(:,:,2),mask);
I2(:,:,3)=klt(I(:,:,3),mask);
figure;
set(gcf,'unit','normalized','position',[0.2,0.5,0.3,0.3]);%设置图形窗大小
set(gcf,'menubar','none'); %删除图形窗工具条
str='量化程度';
set(gcf,'name',[str sprintf('%1.4g\',z)],'numbertitle','off'); %书写图形窗名
imshow(I2);
title('KLT变换压缩后的图像')
imwrite(I2,'klt.jpg');%输出压缩后的图像，文件名为klt.jpg
a=toc;
uicontrol('style', 'text', 'unit', 'normalized', ...
'position', [0 0 0.2 0.1], 'string', num2str(a), ...
'FontSize', 18)
function jpeg_show(I)
% I=imread('football.jpg');
figure;
set(gcf,'unit','normalized','position',[0,0.5,0.3,0.3]);%设置图形窗大小
set(gcf,'menubar','none'); %删除图形窗工具条str='KLT DCT性能比较';
set(gcf,'name',str ,'numbertitle','off'); %书写图形窗名
imshow(I) %显示原始图像
title('原始图像')
%求一个8X8的矩阵的特征矩阵，被调用函数
k=1;
for n1=1:8 %扫描8X8矩阵
for n2=1:8
fz(k)=A(n1,n2);
k=k+1;
end
end%扫描8X8矩阵
k=k-1;
aa=64;
for k=1:64 %将元素变换成一维
X(k,1)=fz(k);
end%将元素变换成一维
for k=1:64 %减去均值后的元素值
fz(k)=fz(k)-jz;
end%减去均值后的元素值
B=zeros(64);
for n1=1:64 %求协方差矩阵
for n2=1:64
B(n1,n2)=fz(n1)*fz(n2);
end
end%求协方差矩阵
%图像性能方差、信息熵及平均梯度。

I=imread('dct48.jpg');
% I=im2double(I);
%%%%%%%%平均梯度
outval = avg_gradient(I)
%%%%%%%%信息熵
S = shannon(I)
%%%%%%%%方差
Q=(std2(I))^2
function outval = avg_gradient(img)
% OUTV AL = A VG_GRADIENT(IMG) 平均梯度
if nargin == 1
img = double(img);
% Get the size of img
[r,c,b] = size(img);
dx = 1;
dy = 1;
for k = 1 : b
band = img(:,:,k);
[dzdx,dzdy] = gradient(band,dx,dy);
s = sqrt((dzdx .^ 2 + dzdy .^2) ./ 2);
g(k) = sum(sum(s)) / ((r - 1) * (c - 1));
end
outval = mean(g);
else
error('Wrong number of input!');
end
function S = shannon(img)
% 信息熵
I=img;
%I=double(I);
[C,R]=size(I); %求图像的规格
Img_size=C*R; %图像像素点的总个数
L=256; %图像的灰度级
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; %求熵值的公式
S=H_img;
end
end。