激光诱导击穿光谱图像有损无损联合压缩方法与流程

本发明属于图像处理技术领域，更进一步涉及一种图像压缩的方法，可用于图像传输。

背景技术：

深空环境等图像传输领域，由于数据传输和存储资源的限制，必须对图像进行压缩处理。目前在图像压缩领域，有一些直接对整幅图进行压缩的方法，如jpeg-ls方法，它是一种新的针对连续色调静态图像的压缩标准(iso-14495-1/itu-t.87)。但在一些特殊应用的需求中，需要保护重要区域的信息免于损失，同时整幅图达到一定的压缩比，因而直接用现有jpeg-ls方法压缩无法满足有损无损同时压缩的需求。

哈尔滨工业大学在其专利申请“一种高光谱图像压缩方法”(专利申请号：cn201110122607.3，公开号：cn102156998a)中公开了一种对高光谱图像进行分级压缩处理的方法。基于jpeg-2000算法实现了对感兴趣空间信息和感兴趣光谱信息无损或者近无损的压缩，而对其它信息进行较大压缩比的压缩。但仍然存在的不足是，此方法需要对图像进行小波变换，实现复杂，压缩速度较低，不能满足对压缩速度要求较高的场合。

现有jpeg-ls标准的核心算法是低复杂度无损图像压缩loco-i结合了哈夫曼编码的低复杂度和上下文建模的高压缩效率，并且压缩速度很高。但现有jpeg-ls方法只能对整幅图进行相同倍数的压缩，有损压缩时整幅图都会有信息损失，在有些需要保证图像特定区域或特定目标信息免于损失的应用中，不能满足要求。

技术实现要素：

本发明的目的在于克服上述现有技术的不足，提出一种激光诱导击穿光谱图像有损无损联合压缩方法，以在保证整幅图像达到设定压缩倍数的前提下，避免图像特定区域或特定目标的信息损失，并且减小压缩复杂度，提高压缩速度。

为实现上述目的，本发明的技术方案包括如下：

(1)提供图像参数：包括待压缩图像、图像宽度w、图像高度h、比特深度dp、需要满足的压缩倍数r和需要无损压缩的区域坐标；

(2)初始化关键变量：

设无损压缩时重构像素与原像素的最大允许误差值为near1，将near1值初始化为0；设有损压缩时重构像素与原像素的最大允许误差值为near2，将near2值初始化为1；

设无损区域的三级量化区间阈值为t1、t2、t3，设有损区域的三级量化区间阈值为t4、t5、t6；

(3)获得当前像素的量化局部梯度值：

(3a)设当前像素正左方的值为b、左上方的值为c、正上方的值为a、右上方的值为d，计算当前像素的3个局部梯度值分别为：d1＝d-b、d2＝b-c、d3＝c-a；

(3b)计算不同区域的量化区间：

根据(2)中设置的t1、t2、t3、near1参数以及它们对应的4个负数值，得到无损压缩的9个量化区间；

根据(2)中设置的t4、t5、t6、near2参数以及它们对应的4个负数值，得到有损压缩的9个量化区间；

(3c)根据(3b)中的量化区间对(3a)中的3个局部梯度值d1、d2、d3进行量化：

若当前像素在无损区域，则根据无损量化区间得到无损的量化局部梯度值q1、q2、q3；

若当前像素在有损区域，则根据有损量化区间得到有损的量化局部梯度值q4、q5、q6；

(4)分别对各区域的量化梯度进行合并，并计算和修正当前像素的预测值，得到无损修正预测值p1和有损修正预测值p2；

(5)根据当前像素值和修正预测值，分别计算得到无损预测误差值err1和有损预测误差值err2；

(6)对预测误差值进行编码：

(6a)计算哥伦布编码参数：

定义无损区域预测误差的出现次数的累积值为n1[q]，无损区域预测误差绝对值的累积值为a1[q]，则无损哥伦布编码参数k1＝|log2(a1[q]/n1[q])|；

定义有损区域预测误差的出现次数的累积值为n2[q]，有损区域预测误差绝对值的累积值为a2[q]，则有损哥伦布编码参数k2＝|log2(a2[q]/n2[q])|；

(6b)在无损区域，根据无损哥伦布编码参数k1，直接在空间域上对无损预测误差值err1进行编码，在有损区域，根据有损哥伦布编码参数k2，直接在空间域上对有损预测误差值err2进行编码，并在编码时限制编码长度，保证压缩性能；

(7)定义累积变量组，并进行参数更新：

定义b1[q]为无损区域修正的预测误差的累积值，并将其与无损区域预测误差绝对值的累积值a1[q]，无损区域预测误差出现次数的累积值n1[q]，无损区域的预测修正值c1[q]统一定义为无损累积变量组；根据无损区域预测误差值err1更新无损累积变量组中的a1[q]、b1[q]、n1[q]，根据更新后b1[q]和n1[q]的值更新c1[q]；

定义b2[q]为有损区域修正的预测误差的累积值，并将其与有损区域预测误差绝对值的累积值a2[q]，有损区域预测误差出现次数的累积值n2[q]，有损区域的预测修正值c2[q]统一定义为有损累积变量组；根据有损区域预测误差值err2更新有损累积变量组中的a2[q]、b2[q]、n2[q]，根据更新后b2[q]和n2[q]的值更新c2[q]；

(8)完成参数更新后，判断是否结束压缩：

若当前像素不是整幅图的最后一个像素，则返回(3)进入下一个像素的压缩；

若当前像素是整幅图的最后一个像素，则判断当前图像压缩倍数是否达到需要满足的压缩倍数r，如果达到，则结束压缩，否则，将有损压缩时重构像素与原像素的最大允许误差值near2加1，无损压缩时重构像素与原像素的最大允许误差值near1值保持为0，返回(2)重新进行关键变量初始化。

本发明与现有技术相比具有以下优点：

第一，本发明由于对无损和有损区域独立进行变量初始化，且在无损和有损区域使用不同的累积变量组并分别进行变量更新，避免了无损和有损压缩的相互干扰，可实现对一幅图像在不同区域进行选择性压缩，即对部分关键区域进行无损压缩，其他部分进行有损压缩，不仅保证了重要信息没有损失，而且使整幅图像达到一定的压缩倍数，减少需要传输的数据量。

第二，本发明由于直接在空间域上对预测误差值进行编码，不进行小波变换等频域变换，因此相比于需要频域变换的有损无损联合图像压缩方法，本发明实现简单，压缩复杂度低，压缩速度快，在低压缩倍数时压缩质量更高。

附图说明

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

图2为本发明中当前像素与相邻像素的位置关系图。

具体实施方式

下面结合附图对本发明的实施进行详细的描述。

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

步骤1，用户根据工程需要提供图像参数。

用户根据工程需要提供的图像参数包括：待压缩图像数据、图像宽度w、图像高度h、图像比特深度dp、需要满足的压缩倍数r和需要无损压缩的区域坐标(x1,y1,x2,y2)，其中x1为无损压缩区域左上角的横坐标，y1为无损压缩区域左上角的纵坐标，x2为无损压缩区域右下角的横坐标，y2为无损压缩区域右下角的纵坐标。

步骤2，判断当前待压缩像素的压缩类型。

将待压缩图像数据逐像素存入开辟的内存空间，设当前待压缩像素的坐标为(x,y)，根据提供的图像参数判断当前待压缩像素是否在需要无损压缩的区域内：

若x1<＝x<＝x2且y1<＝y<＝y2，则当前像素在需要无损压缩的区域内，需要对该像素进行无损压缩；

否则，当前像素不在需要无损压缩的区域内，需要对该像素进行有损压缩。

步骤3，对压缩过程中的关键变量做初始化工作。

(3a)将无损压缩时重构像素与原像素的最大允许误差near1初始化为0，将有损压缩时重构像素与原像素的最大允许误差值near2初始化为1；

(3b)根据near1计算无损区域的三级量化区间阈值t1、t2、t3，以及无损区域预测误差范围range1和无损区域预测误差值所需比特数bpp1：

(3b1)设图像的最大像素值为maxval，则无损区域预测误差范围range1＝maxval+1，无损区域预测误差值所需比特数bpp1＝log2(range1)；

(3b2)计算t1到t3的初始化值：

t1＝clamp(f*1+2+3*near1,near1+1)

t2＝clamp(f*4+3+5*near1,t1)

t3＝clamp(f*17+4+7*near1,t2)；

其中，clamp(i,j)在i大于最大像素值maxval或i小于j时取值为j，否则取值为i；

f为影响因子，在maxval大于128时，f取值为(min(maxval,4095)+128)/256，否则，取值为256/(maxval+1)；*表示乘法运算；

(3c)根据near2计算有损区域的三级量化区间阈值t4、t5、t6，以及有损区域预测误差范围range2和有损区域预测误差值所需比特数bpp2：

(3c1)设有损区域预测误差范围range2＝|(maxval+2*near2)/(2*near2+1)|+1，有损区域预测误差值所需比特数bpp2＝log2(range2)；

(3c2)计算t4到t6的初始化值：

t4＝clamp(f*1+2+3*near2,near2+1)

t5＝clamp(f*4+3+5*near2,t4)

t6＝clamp(f*17+4+7*near2,t5)

其中，clamp(i,j)在i大于最大像素值maxval或i小于j时取值为j，否则取值为i；影响因子f在maxval大于128时取值为(min(maxval,4095)+128)/256，否则取值为256/(maxval+1)。

步骤4，根据当前像素的值及其相邻位置像素的值，计算当前像素的三个局部梯度值d1、d2、d3；

当前像素与相邻像素的位置关系见图2，其中，rx表示当前像素的值，ra、rb、rc和rd表示相邻像素的值；

由所述的相邻像素的值计算得到如下局部梯度值：

第一个局部梯度值：d1＝rd-rb，

第二个局部梯度值：d2＝rb-rc，

第三个局部梯度值：d3＝rc-ra。

步骤5，对当前像素的三个局部梯度值d1、d2、d3进行量化。

(5a)根据步骤3设置的无损区域的三级量化区间阈值t1、t2、t3、无损压缩时重构像素与原像素的最大允许误差near1以及它们对应的负数值-t1、-t2、-t3、-near1，得到无损压缩的9个量化区间，如表1中第一列所示；

分别根据三个局部梯度值d1、d2、d3所在的无损量化区间，得到三个无损的局部梯度量化值q1、q2、q3，其取值范围为-4到4，如表1中第三列所示；

(5b)根据步骤3设置的有损区域的三级量化区间阈值t4、t5、t6、有损压缩时重构像素与原像素的最大允许误差near2以及它们对应的负数值-t4、-t5、-t6、-near2，得到有损压缩的9个量化区间，如表1中第二列所示；

分别根据三个局部梯度值d1、d2、d3所在的有损量化区间，得到三个有损的局部梯度量化值q4、q5、q6，取值范围为-4到4，如表1中第三列所示。

表1量化区间及其对应的量化值

步骤6，分别对各区域的局部梯度量化值进行合并。

(6a)对无损区域局部梯度量化值的合并：

(6a1)判断无损局部梯度量化值q1、q2、q3是否取反，并设置无损反转标志值sign1：

若无损的局部梯度量化值q1、q2、q3的第一个非零元素是负值，则将q1、q2、q3的值取反，并将无损反转标志值sign1设置为-1；

否则，q1、q2、q3保持不变，并将无损反转标志值sign1设置为1，以保证q1、q2、q3的第一个非零元素为正值，取值范围为0到4，有利于局部梯度量化值的合并；

(6a2)按如下条件对q1、q2、q3进行合并，得到无损的合并局部梯度值q：

如果q1和q2都为0，则合并局部梯度值q＝q3+360，其取值范围为360到364；

如果q1为0，而q2不为0，则合并局部梯度值q＝324+(q2-1)*9+(q3+4)，其取值范围为324到359；

如果q1和q2都不为0，则合并局部梯度值q＝(q1-1)*81+(q2+4)*9+(q3+4)，其取值范围为0到323，其中*表示乘法运算；

(6b)对有损区域局部梯度量化值的合并：

(6b1)判断有损局部梯度量化值q4、q5、q6是否取反，并设置有损反转标志值sign2：

若有损的局部梯度量化值q4、q5、q6的第一个非零元素是负值，则将q4、q5、q6的值取反，并将无有损反转标志值sign2设置为-1；

否则，q4、q5、q6保持不变，并将有损反转标志值sign2设置为1，以保证q4、q5、q6的第一个非零元素为正值，取值范围为0到4，有利于局部梯度量化值的合并；

(6b2)按如下条件对q4、q5、q6进行合并，得到有损的合并局部梯度值q，

如果q4和q5都为0，则合并局部梯度值q＝q6+360，其取值范围为360到364；

如果q4为0且q5不为0，则合并局部梯度值q＝324+(q5-1)*9+(q6+4)，其取值范围为324到359；

如果q4和q5都不为0，则合并局部梯度值q＝(q4-1)*81+(q5+4)*9+(q6+4)，其取值范围为0到323。

步骤7，计算得到当前像素的预测值并对其进行修正。

(7a)根据当前像素相邻位置的像素值，使用现有的边缘检测预测器，计算当前像素的预测值px，即根据图2中ra、rb、rc的值来计算px，具体计算公式如下：

(7b)对当前像素的预测值px进行修正：

(7b1)在无损区域，设无损预测修正参数为c1[q]，用c1[q]修正当前像素的预测值px得到无损修正预测值p1，具体修正方法如下：

若步骤(6a1)中无损反转标志值sign1为1，则p1＝px+c1[q]，

若无损反转标志值sign1为-1，则p1＝px-c1[q]；

将p1的值限制在0到maxval范围内，maxval为图像的最大像素值，判断p1的值是否超出修正预测值的上、下限范围：

若p1的值超出修正预测值范围的上限maxval，则将p1的值更改为maxval，

若p1的值超出修正预测值范围的下限0，则将p1的值更改为0；

(7b2)在有损区域，设有损预测修正参数为c2[q]，用c2[q]修正当前像素的预测值px得到有损修正预测值p2，具体修正方法如下：

若步骤(6b1)中有损反转标志值sign2为1，则p2＝px+c2[q]，

若有损反转标志值sign2为-1，则p2＝px-c2[q]；

将p2的值限制在0到maxval范围内，判断p2的值是否超出修正预测值的上、下限范围：

若p2的值超出修正预测值范围的上限maxval，则p2的值更改为maxval，若p2的值超出修正预测值范围的下限0，则p2的值更改为0。

步骤8，计算得到预测误差值。

(8a)计算得到无损预测误差值，并对其进行约束：

(8a1)由当前像素值rx减去无损修正预测值p1得到无损预测误差值err1，并根据步骤(6a1)中无损反转标志值sign1的计算结果决定err1是否取反：若sign1的值为-1，则需要对err1取反，若sign1的值为1，则err1保持不变；

(8a2)根据无损压缩时重构像素与原像素的最大允许误差near1更改无损预测误差值err1：当err1大于0时，err1＝(err1+near1)/(2*near1+1)，否则，err1＝(err1–near1)/(2*near1+1)，其中*表示乘法运算；

(8a3)将无损预测误差值err1取模约束在[-range1/2,range1/2)范围内，其中range1为无损区域预测误差范围，其值在步骤(3b1)中进行了初始化；

(8b)计算得到有损预测误差值，并对其进行约束：

(8b1)由当前像素值rx减去有损修正预测值p2，得到有损预测误差值err2，并根据步骤(6b1)中有损反转标志值sign2的计算结果决定err2是否取反：若sign2的值为-1，则需要对err2取反，若sign2的值为1，则err2保持不变；

(8b2)根据有损压缩时重构像素与原像素的最大允许误差near2更改有损预测误差值err2：当err2大于0时，err2＝(err2+near2)/(2*near2+1)，否则，err2＝(err2–near2)/(2*near2+1)；

(8b3)将有损预测误差值err2取模约束在(-range2/2,range2/2)范围内，其中range2为有损区域预测误差范围，其值在步骤(3c1)中进行了初始化。

步骤9，对预测误差值进行编码。

(9a)在无损区域，对无损预测误差值err1进行编码：

(9a1)计算无损区域哥伦布编码参数：

定义无损区域预测误差的出现次数的累积值为n1[q]，无损区域预测误差绝对值的累积值为a1[q]，则无损哥伦布编码参数k1＝|log2(a1[q]/n1[q])|；

(9a2)根据无损哥伦布编码参数k1，直接在空间域上对无损预测误差值err1进行编码，并在编码时限制编码长度，保证压缩性能，具体编码方法如下：

计算第一中间变量m1：m1等于err1除以2的k1次方；

计算第二中间变量m2：若m1的值大于limit-bpp1-1，将m1的值更改为limit-bpp1-1，第二中间变量m2等于err1-1二进制值的低bpp1位组成的数，否则，第二中间变量m2等于err1二进制值的低k1位组成的数；其中limit为编码长度限定值，bpp1为无损区域预测误差值所需比特数；

由m1、m2组成编码值：将m1个0，一个1，以及m2的二进制值依次排列，形成一串由二进制数0、1组成的值，并将其作为无损预测误差值err1的编码值；

(9b)在有损区域，对有损预测误差值err2进行编码：

(9b1)计算有损区域哥伦布编码参数：

定义有损区域预测误差的出现次数的累积值为n2[q]，有损区域预测误差绝对值的累积值为a2[q]，则有损哥伦布编码参数k2＝|log2(a2[q]/n2[q])|；

(9b2)根据有损哥伦布编码参数k2，直接在空间域上对有损预测误差值err2进行编码，并在编码时限制编码长度，保证压缩性能，具体编码方法如下：

计算第三中间变量m3：m3等于err2除以2的k2次方；

计算第四中间变量m4：若m3的值大于limit-bpp2-1，将m3的值更改为limit-bpp2-1，第四中间变量m4等于err2-1二进制值的低bpp2位组成的数，否则，第四中间变量m4等于err2二进制值的低k2位组成的数，其中bpp2为有损区域预测误差值所需比特数；

由m3、m4组成编码值：将m3个0，一个1，以及m4的二进制值依次排列，形成一串由二进制数0、1组成的值，并将其作为有损预测误差值err2的编码值。

步骤10，定义累积变量组，并对其进行参数更新。

(10a)在无损区域，定义无损累积变量组并更新其值：

(10a1)定义b1[q]为无损区域修正的预测误差的累积值，并将其与无损区域预测误差绝对值的累积值a1[q]、无损区域预测误差出现次数的累积值n1[q]、无损区域的预测修正值c1[q]统一定义为无损累积变量组；

(10a2)根据无损区域预测误差值err1更新无损累积变量组中的a1[q]、b1[q]、n1[q]，根据更新后b1[q]和n1[q]的值更新c1[q]，具体更新如下：

将无损区域预测误差绝对值的累积值a1[q]更新为a1[q]＝a1[q]+abs(err1)，其中abs(.)表示取绝对值；

将无损区域修正的预测误差的累积值b1[q]更新为b1[q]＝b1[q]+err1*(2*near1+1)，其中*表示乘法运算；

将无损区域预测误差出现次数的累积值n1[q]更新为n1[q]＝n1[q]+1；

当更新后的b1[q]小于等于更新后n1[q]的相反数时，将c1[q]更新为c1[q]＝c1[q]-1；当更新后的b1[q]大于0时，将c1[q]更新为c1[q]＝c1[q]+1。

(10b)在有损区域，定义有损累积变量组并更新其值：

(10b1)定义b2[q]为有损区域修正的预测误差的累积值，并将其与有损区域预测误差绝对值的累积值a2[q]，有损区域预测误差出现次数的累积值n2[q]，有损区域的预测修正值c2[q]统一定义为有损累积变量组；

(10b2)根据有损区域预测误差值err2更新有损累积变量组中的a2[q]、b2[q]、n2[q]，根据更新后b2[q]和n2[q]的值更新c2[q]，具体更新如下：

将有损区域预测误差绝对值的累积值a2[q]更新为a2[q]＝a2[q]+abs(err2)，其中abs(.)表示取绝对值；

将有损区域修正的预测误差的累积值b2[q]更新为b2[q]＝b2[q]+err1*(2*near1+1)，其中*表示乘法运算；

将有损区域预测误差出现次数的累积值n2[q]更新为n2[q]＝n2[q]+1；

当更新后的b2[q]小于等于更新后的n2[q]的相反数时，将c2[q]更新为c2[q]＝c2[q]-1；当更新后的b2[q]大于0时，将c2[q]更新为c2[q]＝c2[q]+1。

步骤11，完成参数更新后，判断是否结束压缩。

若当前像素不是整幅图的最后一个像素，则返回步骤4进入下一个像素的压缩；

若当前像素是整幅图的最后一个像素，则判断当前图像压缩倍数是否达到需要满足的压缩倍数r，如果达到，则结束压缩，否则，将有损压缩时重构像素与原像素的最大允许误差值near2加1，无损压缩时重构像素与原像素的最大允许误差值near1值保持为0，返回步骤3重新进行关键变量初始化。

本发明的效果通过以下仿真实验进一步说明：

一.仿真条件

仿真使用的图像为1000张激光诱导击穿光谱图像，图像宽度为64像素、图像高度为64像素、比特深度为12比特；

使用峰值信噪比psnr来评价图像压缩效果，psnr是用来衡量经过压缩处理后图像质量的指标，psnr值越大，图像经过压缩处理后的质量越高，其计算公式如下：

其中mse为仿真图像与恢复图像之间均方误差，n为图像的比特深度，恢复图像为仿真图像经过压缩处理再进行对应的解压操作之后得到的图像。

二.仿真内容

仿真1.用本发明对1000张仿真图像进行压缩，要求对图像中心宽度为16像素且高度为16像素的矩形区域进行无损压缩，且整幅图像依次达到2倍、3倍、4倍、5倍的压缩倍数，完成压缩后，对不同压缩倍数下的压缩码流进行对应的解压操作得到不同压缩倍数时的恢复图像，根据仿真图像和恢复图像计算每幅图的峰值信噪比并对1000幅图的结果取平均值，得到不同压缩倍数时的平均峰值信噪比，具体结果见表1第一列。

仿真2.用现有的基于jpeg_ls标准的图像压缩方法对1000张仿真图像进行压缩，要求整幅图像依次达到2倍、3倍、4倍、5倍的压缩倍数；完成压缩后，对不同压缩倍数下的压缩码流进行对应的解压操作得到不同压缩倍数时的恢复图像，根据仿真图像和恢复图像计算每幅图的峰值信噪比并对1000幅图的结果取平均值，得到不同压缩倍数时的平均峰值信噪比，具体结果见表1第二列。

仿真3.用现有的一种高光谱图像压缩方法对1000张仿真图像进行压缩，要求对图像中心宽度为16像素且高度为16像素的矩形区域进行无损压缩，且整幅图像依次达到2倍、3倍、4倍、5倍的压缩倍数，完成压缩后，对不同压缩倍数下的压缩码流进行对应的解压操作得到不同压缩倍数时的恢复图像，根据仿真图像和恢复图像计算每幅图的峰值信噪比并对1000幅图的结果取平均值，得到不同压缩倍数时的平均峰值信噪比，具体结果见表1第三列。

表1本发明与现有技术的综合比较

从表1可以看出，本发明与现有的jpeg_ls图像压缩方法相比，可以支持部分关键区域无损压缩且其余部分有损压缩的有损无损联合图像压缩，现有jpeg_ls图像压缩方法只能对图像进行整体无损压缩或整体有损压缩。在2倍压缩时，两种方法峰值信噪比基本相同，随着压缩倍数的增大，本发明相比现有jpeg_ls图像压缩方法峰值信噪比有所降低，但差别不大，说明本发明的压缩质量相比jpeg_ls图像压缩方法没有明显下降。

从表1可以看出，本发明与现有的一种高光谱图像压缩方法相比，虽然都可以支持部分关键区域无损压缩且其余部分有损压缩的有损无损联合图像压缩，但在低于4倍压缩时，本发明的峰值信噪比大于一种高光谱图像压缩方法，说明在低压缩倍数时，本发明的压缩质量高于一种高光谱图像压缩方法。