JPEG2000中如何计算失真的？

 

JPEG2000中如何计算失真的？

概述

JPEG2000失真的计算是其EBCOT算法的基础，因此了解如何计算失真才能真正理解EBCOT算法。

首先描述wmse的计算公式：

Wmse  = ( Delta *1<<K_max)^2 * G_b * W_b^2 。

Pass_wmse = wmse/1<<32

Pass_wmse = Pass_wmse * (0.25^miss_msb)

再描述失真的公式：

Delta_D = Pass_wmse*(Ts +Tm)

 

失真的计算牵涉到很多概念，主要包括：

1． 能量计算(G_b)

2． 解码的差异

3． 子带的权重(W_b)

4． Delta的值(Delta) 和K_max

5． Miss_msb和Pass_wmse

6． Ts和Tm

 

能量的计算

首先看一下能量的计算；能量的计算与选择的DWT小波核有关系，这里以5/3小波为例。

看一下5/3小波变换的过程：

正变换公式:

C(2i+1) = P(2i+1) + (1 - P(2i) - P(2i+2))>>1           (1)

C(2i) = P(2i) + (C(2i-1) + C(2i+1) + 2)>>2             (2)

 

反变换公式:

Q(2i) = C(2i) - (C(2i-1) + C(2i+1) + 2)>>2              (3)

Q(2i+1) = C(2i+1) - (1 - Q(2i) - Q(2i+2))>>1           (4)

 

在JPEG2000中对于这种变化叫做提升；5/3变换的提升步骤有两步，分别对应公式1和2，其提升因子分别是-0.5和0.25。

 

同样JPEG2000使用单位脉冲响应来计算能量的变换；这种能量的变换对于高频和低频是不同的。分别在奇数位置和偶数位置设置一个1来模拟高频和低频脉冲响应。下表演示了在整合过程中单位脉冲的响应：这里假设DWT Lelvel为3。

	L1 cof	L1	L2 cof	L2	L3 cof	L3
0
1		0.5		0.25		0.125
2	1	1	0.5	0.5	0.25	0.25
3		0.5		0.75		0.325
4			1	1	0.5	0.5
5				0.75		0.625
6			0.5	0.5	0.75	0.75
7				0.25		0.875
8					1	1
9						0.875
10					0.75	0.75
11						0.625
12					0.5	0.5
13						0.325
14					0.25	0.25
15						0.125

上表是一个在偶数位置模拟一个脉冲的过程，也就是在低频子带中。

Lx表示DWT的第几层

Lx cof表示DWT的整合系数；第一次的整合系数是1，在综合之后得到L1层真正DWT系数；而L1层的系数作为L2层的输入系数，因为有一个放大的过程，所以对应的位置需要乘以，也就是偶数位置，实际上对应低频子带。

由于5/3变换是整数变化，情况和上面有点不同。

 

JPEG2000中认为DWT变换是线性的，实际上从上面的计算可以发现这个事实；因此它认为各层之间的能力是可以线性相加的，从而得到整个过程的能量。

能量是通过系数的平方累加得到。

再来看高频子带中的一个模拟脉冲情况：

	L1 cof	L1	L2 cof	L2	L3 cof	L3
-5
-4
-3
-2		0	-0.125	-0.125
-1	0	-0.125	0
0	0	-0.25	-0.25
1	1	0.75	0
2	0	-0.25	0.75
3	0	-0.125	0
4	0	0	-0.25
5	0		0
6			-0.125
7			0
8
9
10
11
12
13
14
15

 

如果当前的LEVEL是1，而且是低频子带，那么能量系数是：

L_G_b = 0.5*0.5 + 1*1 + 0.5*0.5 = 1.5

如果是高频子带，那么能力系数是：

H_G_b = (-0.125)^2 * 2 + 0.75^2 + (-0.25)^2 = 0.03125+0.5625+0.125=0.71875

 

现在知道如何计算单位冲击的能量了，这些系数可以作为其他系数变化的基数。

 

但DWT变换中某些子带兼有高频和低频部分，可以通过这个公式来计算某个子带的具体能量：

假设子带具有索引x和y，0表示低频，1表示高频，因此x=0，y=0表示LL子带，x=1并且y=1表示HH子带，那么：

G_b = (x? H_G_b: L_G_b) * (y? H_G_b: L_G_b)

分配到子带就是：

子带	G_b
LL	2.25
HL/LH	1.078125
HH	0.5166015625

 

 

不同的分量采用不同的G_b系数：

Y分量：Y_G_b = G_b

U和V分量：G_b = G_b * (0.75^2+0.25^2+0.25^2)

子带的权重

JPEG2000可以为不同的子带分配不同的权重。理论上可以根据CSF来分配权重，但这个值一般都设置为1。

 

Delta和K_max

Delta对于5/3变换实际采用的一个固定的值1/256.

K_max用来表示子带样本范围，不同子带范围是不同的，基本上这样定义：

LL: K_max ＝ bit ＋ B1 ＋1-G

HL/LH: K_max ＝ bit + B2 +1-G

HH: K_max ＝ bit + B3 + 1 – G

G为保护位，一般取1，而B1到B3是最坏的BIBO位扩展，根据计算得到，B1从来不大于2，B2从来不大于3，B3从来不大于4。因此有：

LL: K_max ＝ bit ＋ 2 ＋1-G = bit +2

HL/LH: K_max ＝ bit + 3 +1-G = bit + 3

HH: K_max ＝ bit + 4 + 1 – G= bit + 4

 

Wmse

有来以上各值，我们看一下WMSE是什么值。

这里假设图像的分量精度是8，那么：

对于不同的子带K_max为：

子带	K_max
LL	10
HL/LH	11
HH	12

Wmse = ( Delta *1<<K_max)^2 * G_b * W_b^2

＝ （1/256 ×K_max）^2 *G_b (W_b=1)

结合上面的子带能量分布有：

子带	WMSE
LL	64*G_b=64*2.25 = 144
HL/LH	256* G_b  = 256*1.078125=276
HH	1024*G_b = 1024*0.5166015625=529

Miss_msb

JPEG2000中采用的是位平面编码，其编码是从不是所有像素都为0的位平面开始编码。例如，8位的像素值；如果每个值都小于128，那么开始编码位就是第6（0为基数）位。通常情况在DWT分析的低层，系数一般比较小，他们在高位都位0，因此使用Miss_msb来记录这个事实，从而节省编码。

Miss_msb表示从最高位开始到第一个不全为0的位平面的位平面数量。

Pass_wmse

JPEG2000使用的是位平面编码，并且将每个位平面分为3个过程编码，分别是清除过程、显著过程合量值改进过程。对于第一个编码的位平面只有清除过程，然后从第二个开始编码的位平面开始，都是有：显著过程、量值改进过程合清除过程组成。因此如果有K个位平面需要编码，那么就有：3K-2个编码过程。

Wmse值是对应到每个编码过程就被称为pass_wmse；这个值主要是计算当前过程编码对整个过程的失真改进的情况。

Pass_wmse的计算公式如下：

Pass_wmse = (wmse /2^32)*0.25^miss_msb

由于采用16位整数运算的，所以首先需要将wmse对应到16位整数的变化，因此除以2^16的平方（因为wmse也是一个平方）；然后对于从编码位平面的最高位开始，每前进一个位能量相当于是原来的1/2（因为从数值上来看，例如：0x100到0x080，确实后者是前者的一半），同样需要取平方值。因此有多少个miss_msb就需要缩小次。

这个pass_wmse再编码过程中也是逐渐减小的，每次编码一个位平面完成以后，减少一次。

 

 

Ts和Tm

首先看一下JPEG2000认为如何计算失真的？

D ＝ G_b ∑(y_p[j]-y[j]) ²

其中第一个y是解码后的值，第二个y是原来的样本值；当然这些都是针对DWT系数的。

这里考察从位平面p+1到p时候失真的改进。

Delta_D ＝ G_b×∑[(y_p+1[j]-y[j]) ² - (y_p[j]-y[j]) ² ]

G_b是能量因子，与我们前面的计算值同一个含义。

y_p+1是解码到p＋1位平面时的值，y_p是编码到位平面p时候的值。而y表示原来的DWT系数。

 

如果采用重点重建的算法，那么有：

如果v_p>0： y_p ＝ sign（y） × 2^p×△（v_p+1/2）

如果v_p=0： y_p＝ 0

因此向下编码到位平面p的时候，由样本y产生的失真为：

G_b ×(y_p[j]-y[j]) ²

= G_b × （y-2^p×△（v_p+1/2））² ，其中v_p>0

＝G_b × （y）² ＝（y－2^p×△×v_p）²  ，其中v_p>0

 

上面各式中y的值都是绝对值。

现则假设：v_p为y/2^p×△的小数部分，也就是说：

其：v_p ＝|y|/2^p×△ –v_p

其中v_p表示|y|/2^p×△的整除结果。

 

因此将上面两个结果代入前面的表达式中。这里分为几种情况来计算：

1． v_p为0，v_p+1都为0，这个时候实际上式没有失真的；这中情况不表示

2． v_p为1，表示v_p+1为0，

3． 如果v_p+1为1，而且v_p为0，表示从p+1到p位平面没有改进

4． 如果v_p大于1，表示v_p+1大于0

下面对2和4来求最中的失真改进。

 

对于2有：

(y_p+1[j]-y[j]) ² ＝ ( |y| － 2^p+1△v_p+1 )² = 2^2(p+1)×△×（|y|/2^p+1△ –v_p+1）² =2^2(p+1)×△×v__p+1²

= 2^2p×△×(2v__p+1)²

(y_p[j] - y[j]) ² ＝ (|y| － 2^p△(v_p+1/2) )² = 2^2p×△×（|y|/2^p△ –(v_p+1/2））²= 2^2p×△×(v__p–1/2)²

 

因此有：

Delta_D ＝ G_b×∑[(y_p+1[j]-y[j]) ² - (y_p[j]-y[j]) ² ]

＝ G_b ×2^2p×△×∑[(2v__p+1)²- (v__p–1/2)²]

 

对于4有：

(y_p+1[j]-y[j]) ²＝ ( |y| － 2^p+1△(v_p+1 +1/2))^{2 =} 2^2p×△×(2v__p+1–1/2)²

(y_p[j] - y[j]) ² ＝ (|y| － 2^p△(v_p+1/2) )² = 2^2p×△×（|y|/2^p△ –(v_p+1/2））²= 2^2p×△×(v__p–1/2)²

Delta_D ＝ G_b×∑[(y_p+1[j]-y[j]) ² - (y_p[j]-y[j]) ² ]

= G_b ×2^2p×△×∑[(2v__p+1–1/2)²- (v__p–1/2)²]

 

有一个事实是不可忽略的，也就是：v_p和v_p+1之间的关系。我们知道v_p是|y|/2^p△的小数部分，也就是：

v_p = |y|/2^p△ – [|y|/2^p△]

而v_p+1是|y|/2^P+1△ – [|y|/2^p+1△]。因此2 ×v_p+1是由最第有效位整数位和v_p构成，也就是：

V_p = 2v_p+1 – [2v_p+1],其中[2v_p+1]表示取整。

 

例如：设y＝22，p＝2，p+1＝3，

那么v_p = float（22/4） － （int）22/4 ＝ 0.5

而v_p+1＝ float（22/8） － （int）22/8 ＝ 0.75

因此v_p＝ 2 ×v_p+1－(int)2×v_p+1 ＝ 1.5－1 ＝ 0.5

 

如果样本在位平面p变成显著位，你们v_p ＝1 并且2×v_p+1 >=1（至少应该是1，因为p位至少是p+1位的1/2）。那么失真可表示位：

Delta_D ＝ G_b ×△²×2^2p×∑ Ts（v_p+1）

其中Ts（v_p+1）可表示为：

（2×v_p+1）²－（v_p -1/2）² ＝（2×v_p+1）² －[（2×v_p+1）－1－1/2]²

因为2×v_p+1 >=1，所以其整数部分肯定是1。

对于v_p大于1的过程，一般是量值改进过程，对应的失真计算方法有：

Delta_D ＝ G_b ×△²×2^2p×∑ Tm（v_p+1）

其中Tm（v_p+1）为：

（2v_p+1-1）² － （（2×v_p+1）－x－1/2）²

其中x表示（int）(v_p+1<<1)。如果p位也为1，那么x＝1，否则x＝ 0。

 

JPEG2000中的Ts和Tm

在JPEG2000中通过两个小的查找表来近似计算；这里取5位作为显著过程和清除过程的计算基数。

Ts的计算：

假设将最中的结果转换到1<<16范围内。

设n大于等于0小于32，

       n＝ 32＋n

       v＝ （double）n/(double)32

       s_lut[n] = (long)((v² – (v – 1.5)²)×f_scale＋0.5）

 

Tm的计算。计算Tm的时候使用的是6位来近似整个失真过程的变化。同样将结果转换到1<<16范围内。

设n大于0小于64。

       ^V＝（double）n/（double）32

       如果n大于等于32，那么

       R_lut[n]＝ （v-1.0）²－（v-1.5）²

       如果n小于32，那么：

       R_lut[n]＝ （v-1.0）²－（v-0.5）²

 

这里的计算位5和6根据精度选择。

 

最中的失真表达式

现在已经知道了所有失真计算的变量，可以用来表示最中的失真表达式了。如下：

对于清除和显著过程为：

Delta_D = Pass_wmse *∑Ts

对于量值改建过程为：

Delta_D = Pass_wmse* +∑Tm

 

失真长度斜率

失真的改进和长度的变化比值就是失真长度改进，使用lambda来表示，如下：

Lambda＝ Delta_D/Delta_L

标准的斜率曲线应该式下中心原点内凹的曲线，但上面的Lambda式一个小数值，不便于调试和观察，JPEG2000中使用对数形式来将其转换到0到65536之间的数，具体转换过程如下：

设logscale＝256/ln²

那么log_lambda ＝logscale×（ln ^lambda－ln2³²）＋65536

Log_lambda是不能小于2和大于65535的整数。

 

 

posted on 2007-06-25 14:47 笨笨 阅读(2530) 评论(3)  编辑 收藏 引用 所属分类: 压缩算法