肖孝軍，陳智斌

（1.貴州健康職業(yè)學院，貴州 銅仁 554300；2.昆明理工大學 理學院，昆明 650500）

隨著大量的醫(yī)學成像設(shè)備投入臨床應用，核磁共振（nuclear magnetic resonance，MR）、電子計算機斷層掃描技術(shù)（computed tomography，CT）以及超聲波等數(shù)字醫(yī)學圖像在醫(yī)學臨床診斷中發(fā)揮著越來越關(guān)鍵的作用。首先醫(yī)療斷層掃描對于感興趣的身體結(jié)構(gòu)部位往往輸出16 至64 幅切片圖像，超聲和血管造影等檢查甚至會生成3 min 以上的視頻序列，這必然導致數(shù)字醫(yī)學圖像數(shù)據(jù)量快速增加。其次，大多數(shù)醫(yī)學圖像需要很長的保存周期，這使得數(shù)字醫(yī)學圖像需要更多的存儲空間。此外，在遠程醫(yī)療等應用環(huán)境下，要求在更窄的通信寬帶條件下實現(xiàn)醫(yī)學圖像高保真度的傳輸，如不進行有效壓縮，醫(yī)學圖像將占據(jù)大量的存儲空間，對傳輸寬帶造成很大的壓力。最后，醫(yī)學圖像也屬于自然圖像范疇，其像素之間同樣存在冗余性，具備壓縮的可行性。

目前國內(nèi)外流行的圖像壓縮標準是Joint Photographic Experts Group（JPEG）壓縮標準，但是在有損壓縮條件下，JPEG 為了減少圖像存儲空間，在壓縮過程中往往丟掉一部分無關(guān)緊要的數(shù)據(jù)，而丟失的數(shù)據(jù)無法恢復，導致解碼圖像出現(xiàn)塊效應、馬賽克和Gibbs 等現(xiàn)象。JPEG2000 比較明顯的優(yōu)點是沒有JPEG 壓縮中的塊效應，但是在對小波系數(shù)的高頻部分進行量化時，導致高頻數(shù)據(jù)一定程度的衰減，使得JPEG2000解碼圖像出現(xiàn)模糊失真，主要表現(xiàn)為壓縮偽影、震蕩等壓縮痕跡。

為了提高低比特小波壓縮的解碼圖像質(zhì)量，國內(nèi)外學者提出了多種后處理技術(shù)。Bourdon P 等[1]提出一種基于偏微分方程JPEG2000 解碼迭代算法，此算法能夠提高解碼圖像質(zhì)量，但在高壓縮比、迭代過多情況下平滑過甚，會磨平解碼圖像許多紋理成分。Fang 等[2]在JPEG2000 圖像的空間域和每一個小波子帶提出一種局部自適應Pentive 濾波，并得到一種新壓縮損失計算方法。文獻[3]提出一種邊緣保護小波編碼E-D模型，能夠很好地消除邊緣區(qū)域的壓縮偽影。文獻[4]提出一種原始圖像偽影檢測技術(shù)，并采用空間智能雙邊濾波抑制偽影技術(shù)。本文針對低比特JPEG2000 圖像的模糊失真提出一種基于總變分后處理方案，提高了解碼圖像的質(zhì)量，減輕了邊緣區(qū)域的壓縮偽影。

1 JPEG2000 標準

JPEG2000[5,6]分為編碼端與壓縮端，其編碼器的框圖如圖1(a)所示。首先將離散變換應用于源圖像數(shù)據(jù)。然后，在形成輸出碼流（比特流）之前，對變換系數(shù)進行量化和熵編碼。解碼器與編碼器相反（圖1(b)）。首先對碼流進行熵解碼，反量化和逆離散變換，從而得到重建的圖像數(shù)據(jù)。

在編碼端，JPEG2000 編碼需要對原始圖像進行預處理。首先，如果輸入給編碼端的圖像內(nèi)存較高，此時單一編碼器不能處理整幅圖像，將原始圖像分割為非重疊的片，這些非重疊的片分別進行后續(xù)的編碼操作。其次，在對子塊進行離散小波變換（Forward discrete wavelet transform，F(xiàn)DWT）之前，通過減去相同量進行DC 電平移位。例如減去128，則一幅8 比特的圖像的像素就變?yōu)閰^(qū)間[-127,128]范圍內(nèi)的系數(shù)。最后，執(zhí)行色度空間轉(zhuǎn)換，將RGB 色度空間轉(zhuǎn)換為 YCbCr空間，本文考慮的是8 比特灰度圖像，不考慮色度轉(zhuǎn)換。

在量化之前對輸入數(shù)據(jù)進行正交變換，是為了消除圖像像素之間的統(tǒng)計相關(guān)性。在JPEG 標準中采用的是離散余弦變換，而JPEG2000 采用的是離散小波變換。Mallet 小波分解算法將圖像分解成塔式結(jié)構(gòu)，每一層小波分解形成4 個子帶，即水平和垂直方向的低頻子帶LL、水平和垂直方向的高頻子帶LH、垂直和水平方向的高頻子帶HL、水平和垂直方向的高頻子帶HH。

進行小波變換后，通過降低小波系數(shù)的精度對每一個小波系數(shù)進行量化，當量化步長為1 且系數(shù)是整數(shù)時，JPEG2000 是無損壓縮，反之是有損壓縮，其量化過程是不可逆的，導致解碼圖像質(zhì)量降低。其量化公式如式(1)所示

其中， ab（ i ,j ）表示子帶b 的小波系數(shù)， qb（ i ,j ）表示子帶b 的量化系數(shù)， Δb是子帶b 的量化步長，sign 是符號函數(shù)。在量化過程中，量化步長越大，丟失的數(shù)據(jù)越多，解碼圖像質(zhì)量越低，其計算公式如下：

對于量化步長可以分成兩部分理解，2Rb將小波系數(shù)歸一化到區(qū)間［ -1 /2,1/2］， 2-εb（1 + μb/211）表示真正的量化步長。

JPEG2000 解碼端是壓縮端的逆過程，首先進行熵解碼得到量化系數(shù)，其次進行反量化操作得到反量化系數(shù)，其反量化公式如下：

其中，r 是任意選擇的溢出參數(shù)，0 1r≤ ＜ ，其作用是盡可能減少因量化取整的不可逆性造成的壓縮損失。最后對反量化系數(shù)進行逆小波變換（Inverse discrete wavelet transform，IDWT）得到解碼數(shù)據(jù)。

2 基于總變分的JPEG2000 解碼算法

2.1 總變分解碼模型

核磁共振圖像壓縮（MRI）是醫(yī)學圖像處理中一個基本問題，因其成像掃描時間長，導致圖像像素強度不均勻且空間分辨率不高，高壓縮比情況下重建圖像存在困難。熵編碼是完全可逆的，JPEG 壓縮痕跡主要是因本量化取整的不可逆導致解碼圖像質(zhì)量下降，其退化模型可以理解為

其中，v 是被壓縮系數(shù)，W 是小波變換，u 是相應的原始圖像，η 是小波系數(shù)的壓縮損失。文獻[2]中，作者通過對多幅圖像進行不同的壓縮比實驗發(fā)現(xiàn)，壓縮損失η 近似地服從高斯分布。由于不同的壓縮比，很難正確估計量化噪聲的分布情況，但是通過大量仿真發(fā)現(xiàn)，可以用高斯概率密度函數(shù)去近似估計噪聲分布，如圖2 所示。

模型(4)是一個不適定反問題，無法從v 中保證u 的存在性、穩(wěn)定性和唯一性，因此，有必要對求解方案添加約束使其正則化?？紤]到二次正則化函數(shù)不能處理不連續(xù)圖像的邊緣區(qū)域，Rudin，Oshe 和Fatemi 提出一種總變分（TV）正則化模型[5]，結(jié)合對噪聲進行分析將其修正為等式約束最優(yōu)化問題

σ2是壓縮損失的方差，引入拉格朗日乘子 λ ＞ 0，得到如下無約束的最優(yōu)化解碼（ROF）模型

2.2 總變分離散模型

定義離散圖像空間 X = RN×M，當u ∈ X，v ∈ X ，則內(nèi)積為

同樣定義圖像對偶空間Y = X × X，若u ∈ Y，v ∈ Y，則相應的內(nèi)積定義為

對于u ∈ X，定義 ?： X →Y 為向前差分離散梯度算子， ? u = （ ?xu ,?yu ）'，其中 ?xu ， ?yu表示二維圖像的水平與垂直方向的梯度算子

引理[6,7]存在u X*∈ ，p Y*∈ ，滿足

其中，關(guān)于p 的最大問題是關(guān)于u 最小化問題的對偶形式，A 為線性算子，F(xiàn)*，G*分別是F 和G 的共軛函數(shù)。且上述最優(yōu)解 u*∈X， p*∈Y與式(12)的鞍點問題最優(yōu)解是等價的。

根據(jù)上述引理，則式(11)的對偶形式如下

式(13)中的散度算子 XY →:div 是梯度算子?的伴隨算子， div-=?*，具體表達如下

其中，

綜上，模型(11)對應地鞍點問題如下

2.3 模型求解

許多學者提出不同的算法求解TV去噪問題，例如對偶算法[8,9]、ADMM 算法[10,11]、牛頓方法[12]、bergman[13]、自適應投影算法[14]等。本文利用Chambolle 提出的原對偶算法[15,16]求解式(17)的鞍點問題，原對偶算法具體求解以下鞍點問題

其中， F: Y → ［0 ，+∞）和 G: X→ ［0，+∞）是凸函數(shù)，A 是線性算子，則原對偶算法的具體迭代格式如下

這里 ?F , ?G 分別是函數(shù)F 和G 的次梯度算子，（ I + τ? F ）-1表示鄰近點運算符

表1 基于總變分JPEG2000 圖像解碼算法

對于懲罰參數(shù)λ 的選擇影響消除量化噪聲與邊緣信息的平衡，一種常見的辦法是根據(jù)偏差原理選擇λ 去匹配被壓縮系數(shù)與小波系數(shù)的標準差2σ 。另一方面，希望每一次迭代使得式(5)中的2σ逐漸變小，所以把參數(shù)選擇歸為如下最優(yōu)化問題

綜上所述，基于總變分JPEG2000 圖像解碼算法如表1 所示。

3 數(shù)值實驗

本節(jié)對兩幅MR 圖像（MRI-1、MRI-2）在不同壓縮比條件下進行實驗仿真，此外將TV-JPEG2000 解碼算法與高斯濾波、中值濾波等去噪算法進行比較，用于實驗的原始圖像如圖3 所示。迭代步長τ = 0.15，ρ = 0.25，最大迭代次數(shù)maxiter = 400，終止條件 tol = 10-4。用信噪比(SNR)量化圖像質(zhì)量。其定義如下

其中，u 是原始圖像，*u 是解碼圖像，一般來說，SNR 越高，圖像質(zhì)量越好，與原圖相似度越高。

在CR=10.6 情況下對MRI-1圖像進行仿真實驗，圖4 為TV 函數(shù)曲線，圖5 為對應的SNR 曲線，可以看出隨著迭代次數(shù)的增大，目標函數(shù)不斷減小并逐漸達到穩(wěn)定值，從側(cè)面說明本文算法的收斂性，SNR 值不斷增加，說明本文算法的有效性。

圖6 是不同后處理解碼方法在CR=10.6 情況下MRI-1解碼圖像局部區(qū)域放大比較圖。從圖可以看出，JPEG2000 解碼圖像與原始圖像相比出現(xiàn)明顯的壓縮痕跡，本文算法相對JPEG2000 解碼圖像、高斯濾波解碼圖像、中值濾波解碼圖像一定程度上提高了解碼圖像質(zhì)量，減輕了壓縮痕跡，視覺效果較好。對應的SNR 值可以從表1 查找，本文提出的算法其SNR 值為24.75，相比JPEG2000 解碼圖像提高1.06(dB)，高出中值濾波解碼0.87(dB)，高出高斯濾波解碼圖像24.61(dB)。

圖7 展示在CR=21.39 的情況下，MRI-2 不同方法的解碼圖像以及對應的邊緣檢測，此處用的是Sobel邊緣檢測算子。從圖可以看出，本文算法不僅在圖像邊緣地區(qū)能夠一定程度減輕壓縮痕跡，同時也能夠保護圖像的邊緣區(qū)域細節(jié)和紋理成分。從而說明本文算法能夠減輕JPEG2000 解碼圖像因量化產(chǎn)生的量化噪聲，相應各種算法的SNR 值可以在表2 中查找。

表2 不同算法在不同壓縮比的SNR 值

4 結(jié)論

本文主要研究消除低比特核磁共振圖像因JPEG2000 壓縮造成的壓縮痕跡問題，首先通過對JPEG2000壓縮中量化噪聲進行分析建立總變分ROF 最優(yōu)化解碼模型。其次利用原對偶算法求解提出的ROF 解碼模型得到后處理核磁共振圖像解碼圖像。最后為驗證本文算法在消除壓縮痕跡方面的有效性，對兩幅MR 圖像在多種壓縮比情況下進行解碼，并與傳統(tǒng)的高斯濾波、中值濾波進行比較。數(shù)值實驗結(jié)果表明，總變分JPEG2000 解碼圖像視覺效果較好，SNR 值相對較高，一定程度上消除了JPEG2000 解碼圖像的壓縮痕跡。雖然總變分方法能夠提高JPEG2000 醫(yī)學解碼圖像的質(zhì)量，但隨著迭代次數(shù)的增加會磨平解碼圖像的細節(jié)信息，SNR 值會出現(xiàn)下降的問題，非局部總變分是解決這個問題的可行手段。此外，許多時候超聲和血管造影等醫(yī)學檢查會生成三維多通道醫(yī)學動態(tài)視頻序列，接下來多通道三維總變分JPEG2000 醫(yī)學視頻解碼是進一步研究的重要課題。