CN106952314A - 基于频域方向滤波器的递归预测图像压缩方法 - Google Patents

Abstract

本发明公开了一种基于频域方向滤波器的递归预测图像压缩方法，主要解决现有技术中纹理方向检测准确度低，时间复杂度高，预测系数数目多的问题。其步骤为：1)设计原型方向滤波器；2)根据原型方向滤波器生成一组方向滤波器；3)读入原始灰度图像数据，并进行纹理方向判断；4)扩展原始灰度图像数据；5)对扩展图像数据以块为单位进行读取；6)根据读取的块数据建立马尔可夫预测模型；7)对读取的块数据进行预测压缩；8)判断是否完成所有块数据的预测，如果是，则结束压缩，否则数据读取地址跳变，返回步骤5)。本发明能灵活准确地检测纹理方向，时间复杂度低，预测系数数目少，提升了压缩性能，可用于对各种数字设备的图像压缩编码。

Description

基于频域方向滤波器的递归预测图像压缩方法

技术领域

本发明涉及属于图像处理技术领域，更进一步涉及一种图像压缩方法，可用于对各种数字设备的图像压缩编码。

背景技术

信息化时代人们越来越依靠计算机获取和利用信息，而数字化后的多媒体信息具有海量数据性，与当前硬件技术所能提供的计算机存储资源和网络带宽之间有很大的差距，因此信息的存储和传输造成了很大的困难。图像信息作为计算机上最重要的资源，对其进行有效的压缩处理无疑是必要的，也会给人们带来巨大的好处。图像压缩不但是各种动态图像压缩、传输的基础，而且还是影响其效果好坏的重要因素。

预测编码作为图像压缩最主要的方法，也是最基本的方法，在一些对图片压缩质量要求比较高的领域经常被采用，比如文本图像和医用图像。预测编码一般包括两部分：一部分是去相关，大多采用预测器，利用已编码的邻域像素灰度值来预测待编码的像素灰度值；另一部分是对预测误差进行熵编码。

预测器一般分为如下三类:

(1)全局预测器：对不同图片或对一幅图片中的不同像素点，均采用一组固定的预测系数对其进行预测；

(2)局部预测器：对图像的不同局部采用不同的预测系数对其进行预测；

(3)全自适应预测器：预测系数随图像像素点的变化而变化。

基于预测的主要方法有静止图像压缩标准JPEG-LS中采用的中值边缘检测MED方法，该方法对一幅图像的像素从上到下、从左到右逐像素进行编码，其本质是对当前区域进行一个简单的梯度计算，然后沿着与边界平行的方向进行预测。对于当前像素X的预测，是先根据邻近像素之间的梯度选择不同的变长编码模式，然后对当前像素进行预测，将当前像素与预测像素相减得到预测差值，最后对预测差值进行熵编码。

相比于JPEG-LS，基于上下文的自适应无损图像编码方法CALIC，使用一个更加复杂的基于梯度的自适应预测方法GAP。该方法根据上下文模板计算当前像素的水平梯度和垂直梯度，得到水平梯度和垂直梯度的差值，然后根据差值来判断是否存在水平或者垂直边缘，最后根据水平或垂直边缘的变化程度，适当地选取相邻像素的权重计算当前像素的预测值。

上述JPEG-LS、CALIC虽然都采用了基于上下文的自适应预测算法对像素进行预测，但由于两者采取的均是较为简单的梯度计算，不能很好地跟踪纹理走向，因此压缩精度有待提升。近几年，很多学者对快速帧内预测模式选择进行了研究，并引入了Canny边缘方向直方图方法，更加充分地利用图像的纹理特性。

2011年肖化超和周检发表的文章“基于LMS的图像自适应预测编码”中提出了一种自适应更新权值的预测方法，电子设计工程,2011.19(4):p.109-112。其本质是对预测中采用的邻近像素不断更新权值，来得到较小的预测差值，该方法虽然预测比较准确，但由于没有充分利用图像的纹理特性，时间复杂度过高。

2012年，焦昶哲在其发表的论文“基于纹理相似度判断的帧内预测快速算法研究”中公开了一种基于纹理相似度判断的快速选择帧内预测模式的方法，西安电子科技大学,2012。该方法虽然能避免全搜索模式，大幅度降低帧内预测的时间复杂度，但是由于在各种预测模式下各个像素的预测均采用线性预测，预测系数数目较多，又增加了待压缩的数据量。

发明内容

本发明的目的在于克服上述已有技术的不足，提供一种基于频域方向滤波器的递归预测图像压缩方法，以降低时间复杂度，减少预测系数的数目和预测差值，实现压缩性能的提升。

实现本发明目的的技术方案是：将具有任意方向选择特性的一组方向滤波器与递归马尔可夫预测模型相结合，先判断图像块的纹理方向，然后利用邻近该方向上的若干参考像素建立递归马尔可夫预测模型，采用最小二乘估计得到预测模型中参考像素的权值系数，具体步骤包括如下：

(1)设计原型方向滤波器：

(1a)利用商用软件matlab中的函数freqspace确定二维频域响应空间[f₁,f₂]；

(1b)设定原型方向滤波器的傅立叶频谱支撑域参数K，其中K取值为正实数；

(1c)以支撑域参数K和频域响应空间[f₁,f₂]限定傅立叶频谱支撑范围，生成原型方向滤波器的频域响应函数H_d；

(1d)利用窗函数设计时域二维有限长单位脉冲响应数字滤波器h，其中h的频域响应函数为H_d；

(1e)利用商用软件matlab中的freqz2函数，根据步骤(1d)中的h生成原型方向滤波器H；

(2)利用双立方插值法将步骤(1e)中的原型方向滤波器H逆时针旋转不同的角度θ，得到一组L个方向的方向滤波器，其中，L取值为8～72之间的正整数；

(3)从自然图像集中读入一幅大小为M×M、后缀为.raw格式的原始灰度图像，其中，M表示原始灰度图像的宽度和高度，×表示相乘操作；

(4)对原始灰度图像数据进行纹理方向判断：

(4a)对原始灰度图像数据进行高通滤波，滤除图像的主要低频分量，得到高频图像数据；

(4b)将高频图像数据均匀分成大小为N×N的高频图像数据块，读取这些高频图像数据块，并将读取次数k初始化为1，其中，N表示高频图像数据块的宽度和高度，取值为4或8，×表示相乘操作；

(4c)利用步骤(2)得到的L个方向的方向滤波器对第k个高频图像数据块进行频域滤波，并且计算每个方向滤波后的傅立叶频谱能量和S₀,S₁,…,S_i,…,S_L-1，其中i＝0,1,…L-1，并检测出原始图像数据块的纹理主方向；

(4d)判断是否完成对所有高频图像数据块的纹理方向读取，如果是，则执行步骤(5)，否则，k自增1，返回步骤(4c)；

(5)对步骤(3)中的原始灰度图像数据进行上下各三行，左右各三列的边界扩展，得到大小为(M+6)×(M+6)的扩展图像数据，并将数据读取地址初始化为扩展图像数据地址起始位置；

(6)从数据读取地址开始，对扩展图像数据以大小为(N+3)×(N+3)块为单位进行读取；

(7)根据读取的块数据建立马尔可夫预测模型：

(7a)依照步骤(4)得到的纹理主方向，对块数据中不包括上边三行、左边三列的各个像素点，利用与该纹理方向上邻近的若干参考像素点建立二维递归马尔可夫预测模型；

(7b)利用最小二乘估计的方法，估计递归马尔可夫预测模型中参考像素的权值系数；

(8)块数据的预测

(8a)对块数据中不包括上边三行、左边三列的各个像素，根据步骤(7)中已建立的递归马尔可夫预测模型进行预测，得到预测值矩阵；

(8b)将块数据中不包括上边三行、左边三列的各个像素灰度值，与其在预测值矩阵中对应的预测值进行相减，得到差值矩阵R，并对差值矩阵R进行离散余弦变换DCT得到变换系数矩阵T，然后对变换系数矩阵T进行量化、熵编码，得到压缩码；

(8c)对压缩码进行解码、逆量化、逆DCT变换，得到恢复后的差值矩阵R′，并将R′与差值矩阵R相加得到重建值矩阵，根据重建值矩阵对扩展图像数据块中的相应地址处像素值进行更新；

(9)判断是否完成所有块数据的预测，如果是，则结束压缩，否则，数据读取地址跳变为下一个块数据的起始地址，返回步骤(6)。

本发明与现有技术相比较，具有如下优点：

1.对图像纹理方向的检测更加准确、灵活。

本发明由于对原始图像进行高通滤波，减少了噪声对纹理方向判断的影响，并且充分利用图像傅立叶频谱特性，设计一组具有方向任意性、易于实现、灵活性强等特点的方向滤波器，克服了现有技术中不能很好的跟踪图像纹理走向的缺点，实现对图像纹理方向更加准确的检测。

2.降低时间复杂度。

本发明由于根据一组方向滤波器，预先判断图像块的纹理方向，提高了选择预测模式的速度，降低了时间复杂度。

3.提升了预测性能及压缩性能。

本发明由于对原始灰度图像数据进行上下各三行，左右各三列的边界扩展，得到大小为(M+6)×(M+6)的扩展图像数据，并对扩展图像数据以大小为(N+3)×(N+3)块为单位进行读取，考虑了对块数据的左边三列上边三行的参考像素重建值的读取，克服了现有技术简单地对图像分块导致的结果不精确的问题，减少了预测残差，提升了预测性能；同时由于本发明中将具有任意方向选择特性的方向滤波器与递归马尔可夫预测模型相结合，同一个块数据中的所有像素的预测参照同一组权值系数，减少了预测系数的数目，能够更好的适应图像局部纹理特性，得到更为精准的像素预测值，减少了预测差值，进一步提升了压缩性能。

附图说明

图1是本发明的实现流程图；

图2是本发明中的原型方向滤波器频谱支撑域图；

图3是本发明中设计的一组方向滤波器；

图4是本发明的图像时域方向性与频域方向性关系图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步的详细说明，实施例用于说明本发明，但不用来限制本发明的范围。

参照图1本发明的实现步骤如下：

步骤1，设计原型方向滤波器。

本发明的方法将一组方向滤波器与递归马尔可夫预测模型相结合，预先检测图像块的纹理方向，再利用递归马尔可夫预测模型进行预测，实现图像压缩，因此首先设计一个原型方向滤波器。目前，设计原型方向滤波器的方法主要包括三类：基于频域的设计方法、基于小波变换的设计方法和基于Gabor小波的设计方法，本发明采用基于频域的设计方法，具体步骤包括如下：

1a)利用商用软件matlab中的函数freqspace确定二维频域响应空间[f₁,f₂]，其中f₁、f₂分别表示频域响应空间中的水平频率和垂直频率，f₁、f₂从低到高的变化代表着频域响应空间中从低频到高频的的不同频率点；

1b)设定原型方向滤波器的傅立叶频谱支撑域参数K，通过K的设定，确定原型方向滤波器的傅立叶频谱支撑域的张开角度α的大小，α＝2arctan(1/K)，其中K取值为正实数；

1c)由于图像的各方向性分量在其傅立叶频谱图中是按照楔形来分布的，因此以支撑域参数K和频域响应空间[f₁,f₂]限定傅立叶频谱支撑范围，生成原型方向滤波器的频域响应函数H_d，该频域响应函数H_d的取值按如下条件确定：

如果((Kf₁-f₂＜0)&&(K×f₁+f₂＜0))||((K×f₁-f₂＞0)&&(K×f₁+f₂＞0))成立，则H_d取值为0，否则，H_d取值为1，其中，&&代表逻辑与操作，||代表逻辑或操作，K取值为正实数；

1d)利用窗函数设计频域响应函数为H_d的时域二维有限长单位脉冲响应数字滤波器h,其中窗函数采用阻带衰减较大、过渡带带宽相对较窄、长度为64的汉宁窗，使傅立叶频谱的主瓣加宽并降低，旁瓣显著较少，以减少频谱泄漏现象；

1e)利用商用软件matlab中的freqz2函数，根据数字滤波器h生成原型方向滤波器H，原型方向滤波器的长度越大，频域滤波的精度越高，对于4×4的块数据，原型方向滤波器的长度取64即可；

原型方向滤波器H的傅里叶频谱支撑域如图2所示，其中ω_x轴为水平频率坐标轴，ω_y轴为竖直频率坐标轴，矩形框与ω_x轴正半轴的交点为(0,π)、与ω_x轴负半轴的交点为(0,-π)，矩形框与ω_y轴的正半轴的交点为(π,0)、与ω_y轴负半轴的交点为(-π,0)，矩形框内夹角为α的区域为原型方向滤波器H的傅里叶频谱支撑域，其中，α＝2arctan(1/K)。

步骤2，生成一组方向滤波器。

参照支撑域参数K，利用双立方插值法将步骤1中生成的原型方向滤波器H依次旋转θ≈i×α的角度，生成一组方向滤波器，其中α＝2arctan(1/K)，α≈180^°/L，i取值为0～L-1范围内的整数，L为8～72之间的正整数；

L的取值不同，则得到不同的一组方向滤波器，以L＝16为例，一组方向滤波器如图3所示，图3中由左到右、由上到下共有16个方向滤波器，分别为Hrotate0、Hrotate1、Hrotate2、Hrotate3、Hrotate4、Hrotate5、Hrotate6、Hrotate7、Hrotate8、Hrotate9、Hrotate10、Hrotate11、Hrotate12、Hrotate13、Hrotate14和Hrotate15。

通过上述步骤1和步骤2，即可完成一组方向滤波器的设计，以下将开始通过步骤3、步骤4对灰度图像进行纹理方向进行检测：

步骤3，读入原始灰度图像数据。

从自然图像集中读入一幅大小为M×M、后缀为.raw格式的原始灰度图像，其中，M表示原始灰度图像的宽度和高度，×表示相乘操作；

步骤4，对原始灰度图像数据进行纹理方向判断。

4a)利用高斯高通滤波器对原始灰度图像数据进行高通滤波，滤除图像的主要低频分量，得到高频图像数据，以减少噪声对纹理方向判断的影响，其中高斯高通滤波器的传输函数为：其中ω₁＝2πf₁，ω₂＝2πf₂，D₀是截止频率，对于不同的灰度图像，D₀取不同的经验值，D(ω₁,ω₂)为距傅立叶变换中心点的距离；

4b)将高频图像数据均匀分成大小为N×N的高频图像数据块，读取这些高频图像数据块，并将读取次数k初始化为1，其中，N表示高频图像数据块的宽度和高度，取值为4或8，×表示相乘操作；

4c)利用步骤2得到的L个方向的方向滤波器对第k个高频图像数据块按如下步骤进行频域滤波：

4c1)对第k个高频图像数据块进行傅立叶变换，得到傅立叶变换矩阵；

4c2)根据傅立叶变换矩阵和L个方向的方向滤波器，得到L个系数矩阵PF_i，其中i＝0,1,…L-1；

4c3)分别将每个系数矩阵PF_i的各个元素进行相加，得到每个方向滤波后的傅立叶频谱能量和S₀,S₁,…,S_i,…,S_L-1；

4d)检测原始图像数据块的纹理主方向：

(4d1)计算每个方向滤波后的傅立叶频谱能量和S₀,S₁,…,S_i,…,S_L-1中的最大值S_max：

S_max＝max{S₀,S₁,…,S_i,…,S_L-1}

(4d2)根据S_max确定图像频域方向性角度即S_max对应的方向滤波器的方向角度为图像频域方向性角度

(4d3)根据图像的空域方向性与其傅立叶频谱方向性存在的正交关系，得到高频图像数据块的纹理主方向的角度为进而得到原始灰度图像数据块的纹理方向的角度为

以与水平方向正半轴逆时针60度的单纹理图像为例，图像时域方向性与频域方向性关系图如图4所示，其中：

图4(a)为与水平方向正半轴成逆时针60度夹角的单纹理图像；

图4(b)为利用步骤2得到的L个方向的方向滤波器对单纹理图像进行滤波后，每个方向滤波器的傅立叶频谱能量和S₀,S₁,…,S_i,…,S_L-1的直方图，其中直方图的横坐标为1～L之间的正整数，分别对应每个方向滤波器，直方图的纵坐标代表每个方向滤波器滤波后的归一化频谱能量的值；

图4(c)为单纹理图像的傅立叶频谱；

图4(d)为图4(b)中最大的归一化频谱能量对应的方向滤波器；

将图4(a)与图4(c)进行对比，可知图像时域方向性与频域方向性成正交关系，将图4(c)与图4(d)进行对比，可证明设计的一组方向滤波器的准确性；

4e)判断是否完成对所有高频图像数据块的纹理方向读取，如果是，则执行步骤5，否则，k自增1，返回步骤4c)；

通过上述步骤3和步骤4，即可完成对原始灰度图像数据进行纹理方向判断，以下步骤将开始对灰度图像进行预测压缩：

步骤5，扩展原始灰度图像数据。

对步骤3中的原始灰度图像数据进行上下各三行，左右各三列的边界扩展，得到大小为(M+6)×(M+6)的扩展图像数据，并将数据读取地址初始化为扩展图像数据的地址起始位置，其中上下各三行，左右各三列的边界扩展，是指以原始灰度图像的边缘像素点为中心，对扩展图像数据中的上下各三行，左右各三列像素点进行镜像赋值；

步骤6，从数据读取地址开始，对扩展图像数据以大小为(N+3)×(N+3)块为单位进行读取。

(6a)初始化一个(N+3)×(N+3)二维矩阵来存储块数据，并将读取次数i初始化为0，由于对不包括上边三行、左边三列的块数据中的像素的预测，需要参考邻近的及已完成预测的块的若干像素重建值，因此在该步骤中，多分配出3×3的矩阵单元来存储邻近的及已完成预测的块的若干像素重建值；

(6b)从数据读取地址开始，读取N+3个数据，并将这些数据赋值给二维矩阵的第一行；

(6c)将数据读取地址增加M+6，并将i自增1；

(6d)判断i是否等于N+3，如果是，则完成块数据的读取，否则，返回步骤(6b)。

步骤7，根据读取的块数据建立马尔可夫预测模型。

7a)依照步骤4得到的纹理主方向，对块数据中不包括上边三行、左边三列的各个像素点，利用与该纹理方向上邻近的若干参考像素点建立二维递归马尔可夫预测模型如下：

u(i,j)＝ρ₁u₁(i,j)+ρ₂u₂(i,j)+ρ₃u₃(i,j)+ρ₄u₄(i,j)+e(i,j),i、j＝0,1,2…N

其中u(i,j)表示坐标为(i,j)的像素点的预测值，u₁(i,j),u₂(i,j),u₃(i,j),u₄(i,j)分别表示四个参考像素的灰度值，ρ₁,ρ₂,ρ₃,ρ₄分别为四个参考像素的权值系数，e为白噪声；

根据设计的一组方向滤波器的方向数目L的取值不同，可选取不同位置的四个参考像素，以L＝8时为例，八种纹理方向选取的四个参考像素的灰度值如下：

在水平0度方向：

u₁(i,j)＝u(i-1,j)

u₂(i,j)＝u(i-2,j)

u₃(i,j)＝u(i-3,j)

u₄(i,j)＝u(i-2,j-1)

在逆时针与x轴正半轴成22.5度夹角的方向：

u₁(i,j)＝u(i,j-1)

u₂(i,j)＝u(i+1,j-1)

u₃(i,j)＝u(i+2,j-1)

u₄(i,j)＝u(i+2,j-2)

在逆时针与x轴正半轴成45度夹角的方向：

u₁(i,j)＝u(i,j-1)

u₂(i,j)＝u(i+1,j-1)

u₃(i,j)＝u(i+1,j-2)

u₄(i,j)＝u(i+2,j-2)

在逆时针与x轴正半轴成67.5度夹角的方向：

u₁(i,j)＝u(i,j-1)

u₂(i,j)＝u(i+1,j-1)

u₃(i,j)＝u(i,j-2)

u₄(i,j)＝u(i+1,j-2)

在逆时针与x轴正半轴成90度夹角的方向：

u₁(i,j)＝u(i,j-1)

u₂(i,j)＝u(i+1,j-1)

u₃(i,j)＝u(i+1,j-2)

u₄(i,j)＝u(i+2,j-2)

在逆时针与x轴正半轴成112.5度夹角的方向：

u₁(i,j)＝u(i,j-1)

u₂(i,j)＝u(i-1,j-1)

u₃(i,j)＝u(i-1,j-2)

u₄(i,j)＝u(i,j-2)

在逆时针与x轴正半轴成135度夹角的方向：

u₁(i,j)＝u(i-1,j)

u₂(i,j)＝u(i-1,j-1)

u₃(i,j)＝u(i,j-1)

u₄(i,j)＝u(i-2,j-2)

在逆时针与x轴正半轴成157.5度夹角的方向：

u₁(i,j)＝u(i-1,j)

u₂(i,j)＝u(i-2,j-1)

u₃(i,j)＝u(i-1,j-1)

u₄(i,j)＝u(i,j-1)

在L＞8时，针对L种纹理方向分别采用二次插值法，计算在纹理方向上由近及远的四个同向像素灰度值，来作为四个参考像素的灰度值；

7b)利用最小二乘估计的方法，对递归马尔可夫预测模型中参考像素的权值系数按如下步骤进行估计：

7b1)将块数据中不包括上边三行、左边三列的各个像素的预测误差平方和E(i,j)表示如下：

E(i,j)＝(u(i,j)-ρ₁u₁(i,j)-ρ₂u₂(i,j)-ρ₃u₃(i,j)-ρ₄u₄(i,j))²

其中u(i,j)表示坐标为(i,j)的像素点的预测值，u₁(i,j),u₂(i,j),u₃(i,j),u₄(i,j)分别表示四个参考像素的灰度值，ρ₁,ρ₂,ρ₃,ρ₄分别为四个参考像素的权值系数，i、j＝1,2,…N，其中N表示块数据中不包括上边三行、左边三列像素的区域的行数及列数，取值为4或8；

7b2)将块数据中不包括上边三行、左边三列的各个像素的预测误差平方和E(i,j)进行相加，得到块数据的预测误差平方和E_total：

其中表示在i＝1,2,…N、j＝1,2,…N范围内对E(i,j)累加；

7b3)根据块数据的预测误差平方和E_total建立方程组：

其中表示将E_total对第一个参考像素的权值系数ρ₁求偏微分，表示将E_total对第二个参考像素的权值系数ρ₂求偏微分表示将E_total对第三个参考像素的权值系数ρ₃求偏微分，表示将E_total对第四个参考像素的权值系数ρ₄求偏微分；

7b4)对步骤(7b3)中的方程组进行求解，得到参考像素权值系数向量进而得到四个参考像素的权值系数ρ₁,ρ₂,ρ₃,ρ₄。

步骤8，对原始图像块数据进行预测压缩。

8a)对块数据中不包括上边三行、左边三列的各个像素，根据步骤7中已建立的递归马尔可夫预测模型进行预测，得到预测值矩阵；

8b)将块数据中不包括上边三行、左边三列的各个像素灰度值，与其在预测值矩阵中对应的预测值进行相减，得到差值矩阵R，并对差值矩阵R进行离散余弦变换DCT得到变换系数矩阵T，然后对变换系数矩阵T进行量化、熵编码，得到压缩码，其中如果将量化步骤省略，可以得到无损压缩的压缩码；

8c)对压缩码进行解码、逆量化、逆DCT变换，得到恢复后的差值矩阵R′，并将R′与差值矩阵R相加得到重建值矩阵，根据重建值矩阵对扩展图像数据块中的相应地址处像素值进行更新。

步骤9，判断是否完成所有块数据的预测，如果是，则结束压缩，否则，数据读取地址跳变为下一个块数据的起始地址，返回步骤(6)。

以上描述仅是本发明的一个具体实例，不构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解本发明的内容和原理后，都可能在不背离本发明原理、结构的情况下，进行形式和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。

Claims (8)

1.基于频域方向滤波器的递归预测图像压缩方法，包括：

(1)设计原型方向滤波器：

(1a)利用商用软件matlab中的函数freqspace确定二维频域响应空间[f₁,f₂]；

(1b)设定原型方向滤波器的傅立叶频谱支撑域参数K，其中K取值为正实数；

(1c)以支撑域参数K和频域响应空间[f₁,f₂]限定傅立叶频谱支撑范围，生成原型方向滤波器的频域响应函数H_d；

(1d)利用窗函数设计时域二维有限长单位脉冲响应数字滤波器h，其中h的频域响应函数为H_d；

(1e)利用商用软件matlab中的freqz2函数，根据步骤(1d)中的h生成原型方向滤波器H；

(2)利用双立方插值法将步骤(1e)中的原型方向滤波器H逆时针旋转不同的角度θ，得到一组L个方向的方向滤波器，其中，L取值为8～72之间的正整数；

(3)从自然图像集中读入一幅大小为M×M、后缀为.raw格式的原始灰度图像，其中，M表示原始灰度图像的宽度和高度，×表示相乘操作；

(4)对原始灰度图像数据进行纹理方向判断：

(4a)对原始灰度图像数据进行高通滤波，滤除图像的主要低频分量，得到高频图像数据；

(4b)将高频图像数据均匀分成大小为N×N的高频图像数据块，读取这些高频图像数据块，并将读取次数k初始化为1，其中，N表示高频图像数据块的宽度和高度，取值为4或8，×表示相乘操作；

(4c)利用步骤(2)得到的L个方向的方向滤波器对第k个高频图像数据块进行频域滤波，并且计算每个方向滤波后的傅立叶频谱能量和S₀,S₁,…,S_i,…,S_L-1，其中i＝0,1,…L-1，并检测出原始图像数据块的纹理主方向；

(4d)判断是否完成对所有高频图像数据块的纹理方向读取，如果是，则执行步骤(5)，否则，k自增1，返回步骤(4c)；

(5)对步骤(3)中的原始灰度图像数据进行上下各三行，左右各三列的边界扩展，得到大小为(M+6)×(M+6)的扩展图像数据，并将数据读取地址初始化为扩展图像数据地址起始位置；

(6)从数据读取地址开始，对扩展图像数据以大小为(N+3)×(N+3)块为单位进行读取；

(7)根据读取的块数据建立马尔可夫预测模型：

(7a)依照步骤(4)得到的纹理主方向，对块数据中不包括上边三行、左边三列的各个像素点，利用与该纹理方向上邻近的若干参考像素点建立二维递归马尔可夫预测模型；

(7b)利用最小二乘估计的方法，估计递归马尔可夫预测模型中参考像素的权值系数；

(8)块数据的预测

(8a)对块数据中不包括上边三行、左边三列的各个像素，根据步骤(7)中已建立的递归马尔可夫预测模型进行预测，得到预测值矩阵；

(8b)将块数据中不包括上边三行、左边三列的各个像素灰度值，与其在预测值矩阵中对应的预测值进行相减，得到差值矩阵R，并对差值矩阵R进行离散余弦变换DCT得到变换系数矩阵T，然后对变换系数矩阵T进行量化、熵编码，得到压缩码；

(8c)对压缩码进行解码、逆量化、逆DCT变换，得到恢复后的差值矩阵R′，并将R′与差值矩阵R相加得到重建值矩阵，根据重建值矩阵对扩展图像数据块中的相应地址处像素值进行更新；

(9)判断是否完成所有块数据的预测，如果是，则结束压缩，否则，数据读取地址跳变为下一个块数据的起始地址，返回步骤(6)。

2.根据权利要求1所述的方法，其中步骤(1c)中的频域响应函数H_d，其取值按如下条件确定：

如果((Kf₁-f₂＜0)&&(K×f₁+f₂＜0))||((K×f₁-f₂＞0)&&(K×f₁+f₂＞0))成立，则H_d取值为0，否则，H_d取值为1，其中，&&代表逻辑与操作，||代表逻辑或操作。

3.根据权利要求1所述的方法，其中步骤(1d)的窗函数，是指采用长度为64的汉宁窗。

4.根据权利要求1所述的方法，其中步骤(2)中逆时针旋转不同的角度θ，是指参照傅立叶频谱支撑域参数K，依次旋转θ≈i×α的角度，其中α＝2arctan(1/K)，i取值为0～L-1范围内的整数，L为8～72之间的正整数。

5.根据权利要求1所述的方法，其中步骤(4c)中检测出原始灰度图像数据块的纹理主方向，按如下步骤进行

(4c1)计算每个方向滤波后的傅立叶频谱能量和S₀,S₁,…,S_i,…,S_L-1中的最大值S_max：

S_max＝max{S₀,S₁,…,S_i,…,S_L-1}

(4c2)根据S_max确定图像频域方向性角度

(4c3)根据图像的空域方向性与其傅立叶频谱方向性存在的正交关系，得到高频图像数据块的纹理主方向的角度为进而得到原始灰度图像数据块的纹理方向的角度为

6.根据权利要求1所述的方法，其中步骤(6)中对扩展图像数据以大小为(N+3)×(N+3)块为单位进行读取，按如下步骤进行：

(6a)初始化一个(N+3)×(N+3)二维矩阵来存储块数据，并将读取次数i初始化为1；

(6b)从数据读取地址开始，读取N+3个数据，并将其赋值给二维矩阵的第一行；

(6c)将数据读取地址增加M+6，并将i自增1；

(6d)判断i是否等于N+3，如果是，则完成块数据的读取，否则，返回步骤(6b)。

7.根据权利要求1所述的方法，其中步骤(7a)中的二维马尔可夫预测模型，表示如下：

u(i,j)＝ρ₁u₁(i,j)+ρ₂u₂(i,j)+ρ₃u₃(i,j)+ρ₄u₄(i,j)+e(i,j),i、j＝0,1,2…N

其中u(i,j)表示坐标为(i,j)的像素点的预测值，u₁(i,j),u₂(i,j),u₃(i,j),u₄(i,j)分别表示四个参考像素的灰度值，ρ₁,ρ₂,ρ₃,ρ₄分别为四个参考像素的权值系数，e(i,j)为白噪声。

8.根据权利要求1所述的方法，其中步骤(7b)中的估计递归马尔可夫预测模型中参考像素的权值系数，按如下步骤进行：

(7b1)将块数据中不包括上边三行、左边三列的各个像素的预测误差平方和E(i,j)表示如下：

E(i,j)＝(u(i,j)-ρ₁u₁(i,j)-ρ₂u₂(i,j)-ρ₃u₃(i,j)-ρ₄u₄(i,j))²

其中u(i,j)表示坐标为(i,j)的像素点的预测值，u₁(i,j),u₂(i,j),u₃(i,j),u₄(i,j)分别表示四个参考像素的灰度值，ρ₁,ρ₂,ρ₃,ρ₄分别为四个参考像素的权值系数，i、j＝1,2,…N，其中N表示块数据中不包括上边三行、左边三列像素的区域的行数及列数，取值为4或8；

(7b2)将块数据中不包括上边三行、左边三列的各个像素的预测误差平方和E(i,j)进行相加，得到块数据的预测误差平方和E_total：

$E_{total}=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}E(i,j)$

其中表示在i＝1,2,…N、j＝1,2,…N范围内对E(i,j)累加；

(7b3)根据块数据的预测误差平方和E_total建立方程组：

$\left\{\begin{array}{c}\∂E_{total}/\∂{\mathrm{\ρ}}_1=0\\ \∂E_{total}/\∂{\mathrm{\ρ}}_2=0\\ \∂E_{total}/\∂{\mathrm{\ρ}}_3=0\\ \∂E_{total}/\∂{\mathrm{\ρ}}_4=0\end{array}\right.$

其中表示将E_total对第一个参考像素的权值系数ρ₁求偏微分，表示将E_total对第二个参考像素的权值系数ρ₂求偏微分表示将E_total对第三个参考像素的权值系数ρ₃求偏微分，表示将E_total对第四个参考像素的权值系数ρ₄求偏微分；

(7b4)对步骤(7b3)中的方程组进行求解，得到参考像素权值系数向量进而得到四个参考像素的权值系数ρ₁,ρ₂,ρ₃,ρ₄。