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 笨笨 閱讀(2598) 評論(3)  編輯 收藏 引用 所屬分類: 壓縮算法