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 ∑(yp[j]-y[j]) 2
其中第一個y是解碼后的值，第二個y是原來的樣本值；當然這些都是針對DWT系數的。
這里考察從位平面p+1到p時候失真的改進。
Delta_D ＝ G_b×∑[(yp+1[j]-y[j]) 2 - (yp[j]-y[j]) 2 ]
G_b是能量因子，與我們前面的計算值同一個含義。
yp+1是解碼到p＋1位平面時的值，yp是編碼到位平面p時候的值。而y表示原來的DWT系數。
如果采用重點重建的算法，那么有：
如果vp>0： yp ＝ sign（y） × 2p×△（vp+1/2）
如果vp=0： yp＝ 0
因此向下編碼到位平面p的時候，由樣本y產生的失真為：
G_b ×(yp[j]-y[j]) 2
= G_b × （y-2p×△（vp+1/2））2 ，其中vp>0
＝G_b × （y）2 ＝（y－2p×△×vp）2 ，其中vp>0
上面各式中y的值都是絕對值。
現則假設：v_p為y/2p×△的小數部分，也就是說：
其：v_p ＝|y|/2p×△ –vp
其中vp表示|y|/2p×△的整除結果。
因此將上面兩個結果代入前面的表達式中。這里分為幾種情況來計算：
1． vp為0，vp+1都為0，這個時候實際上式沒有失真的；這中情況不表示
2． vp為1，表示vp+1為0，
3． 如果vp+1為1，而且vp為0，表示從p+1到p位平面沒有改進
4． 如果vp大于1，表示vp+1大于0
下面對2和4來求最中的失真改進。
對于2有：
(yp+1[j]-y[j]) 2 ＝ ( |y| － 2p+1△vp+1 )2 = 22(p+1)×△×（|y|/2p+1△ –vp+1）2 =22(p+1)×△×v_p+12
= 22p×△×(2v_p+1)2
(yp[j] - y[j]) 2 ＝ (|y| － 2p△(vp+1/2) )2 = 22p×△×（|y|/2p△ –(vp+1/2））2= 22p×△×(v_p–1/2)2
因此有：
Delta_D ＝ G_b×∑[(yp+1[j]-y[j]) 2 - (yp[j]-y[j]) 2 ]
＝ G_b ×22p×△×∑[(2v_p+1)2- (v_p–1/2)2]
對于4有：
(yp+1[j]-y[j]) 2＝ ( |y| － 2p+1△(vp+1 +1/2))2 = 22p×△×(2v_p+1–1/2)2
(yp[j] - y[j]) 2 ＝ (|y| － 2p△(vp+1/2) )2 = 22p×△×（|y|/2p△ –(vp+1/2））2= 22p×△×(v_p–1/2)2
Delta_D ＝ G_b×∑[(yp+1[j]-y[j]) 2 - (yp[j]-y[j]) 2 ]
= G_b ×22p×△×∑[(2v_p+1–1/2)2- (v_p–1/2)2]
有一個事實是不可忽略的，也就是：v_p和v_p+1之間的關系。我們知道v_p是|y|/2p△的小數部分，也就是：
v_p = |y|/2p△ – [|y|/2p△]
而v_p+1是|y|/2P+1△ – [|y|/2p+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變成顯著位，你們vp ＝1 并且2×v_p+1 >=1（至少應該是1，因為p位至少是p+1位的1/2）。那么失真可表示位：
Delta_D ＝ G_b ×△2×22p×∑ Ts（v_p+1）
其中Ts（v_p+1）可表示為：
（2×v_p+1）2－（v_p -1/2）2 ＝（2×v_p+1）2 －[（2×v_p+1）－1－1/2]2
因為2×v_p+1 >=1，所以其整數部分肯定是1。
對于vp大于1的過程，一般是量值改進過程，對應的失真計算方法有：
Delta_D ＝ G_b ×△2×22p×∑ Tm（v_p+1）
其中Tm（v_p+1）為：
（2v_p+1-1）2 － （（2×v_p+1）－x－1/2）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)((v2 – (v – 1.5)2)×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）2－（v-1.5）2
如果n小于32，那么：
R_lut[n]＝ （v-1.0）2－（v-0.5）2
這里的計算位5和6根據精度選擇。
最中的失真表達式
現在已經知道了所有失真計算的變量，可以用來表示最中的失真表達式了。如下：
對于清除和顯著過程為：
Delta_D = Pass_wmse *∑Ts
對于量值改建過程為：
Delta_D = Pass_wmse* +∑Tm
失真長度斜率
失真的改進和長度的變化比值就是失真長度改進，使用lambda來表示，如下：
Lambda＝ Delta_D/Delta_L
標準的斜率曲線應該式下中心原點內凹的曲線，但上面的Lambda式一個小數值，不便于調試和觀察，JPEG2000中使用對數形式來將其轉換到0到65536之間的數，具體轉換過程如下：
設logscale＝256/ln2
那么log_lambda ＝logscale×（ln lambda－ln232）＋65536
Log_lambda是不能小于2和大于65535的整數。