基于梯度場(chǎng)變換的MRI圖像增強(qiáng)方法

吳夢(mèng)飛， 薛旭成， 蘭太吉， 徐鑫偉

1. 中國(guó)科學(xué)院長(zhǎng)春光學(xué)精密機(jī)械與物理研究所， 吉林 長(zhǎng)春 130033；2. 中國(guó)科學(xué)院大學(xué)， 北京 100049

1 引言

MRI的全稱(chēng)是磁共振成像(magnetic resonance imaging)，在日常醫(yī)學(xué)診斷中，MRI常用于腫瘤、炎癥、退行性病變等疾病的檢查。然而，由于磁共振機(jī)器在工作時(shí)會(huì)產(chǎn)生非常強(qiáng)大的磁場(chǎng)，其所成像的MRI圖像存在較多的系統(tǒng)噪聲，而且部分MRI圖像整體昏暗，細(xì)節(jié)特征表現(xiàn)不清。這無(wú)疑給后續(xù)的診斷工作帶來(lái)了很大的困難。針對(duì)于這一問(wèn)題，目前有學(xué)者提出各種方法用于MRI增強(qiáng)，但這些方法在增強(qiáng)的同時(shí)也會(huì)對(duì)系統(tǒng)噪聲產(chǎn)生過(guò)度增強(qiáng)的效果[1-4]。本文提出一種基于梯度場(chǎng)變換的圖像增強(qiáng)方法，該方法在對(duì)MRI增強(qiáng)的同時(shí)可以有效地避免圖像背景噪聲的增強(qiáng)。實(shí)驗(yàn)結(jié)果表明，本文提出的方法增強(qiáng)效果優(yōu)于直方圖均衡化等傳統(tǒng)的方法。

直方圖均衡化是一種常用的圖像增強(qiáng)算法，它能夠?qū)D像區(qū)域中的灰度范圍進(jìn)行拉伸，使其以相同的概率分布在整個(gè)直方圖上，在醫(yī)學(xué)影像領(lǐng)域得到了廣泛應(yīng)用。對(duì)于對(duì)比度較低的低照度圖像，直方圖均衡化技術(shù)能夠有效地提升其對(duì)比度，產(chǎn)生較好的視覺(jué)效果。另外，因?yàn)橹狈綀D均衡化在理論上等同于熵最大化，它是熵最大化理論的直接體現(xiàn)，這意味著圖像中的特征信息會(huì)以相同的概率表現(xiàn)出來(lái)[5]。但直方圖均衡化也有其局限性，即容易對(duì)圖像產(chǎn)生過(guò)度增強(qiáng)的效應(yīng)。經(jīng)過(guò)直方圖均衡化處理過(guò)的圖像，其部分細(xì)節(jié)特征會(huì)產(chǎn)生過(guò)亮的現(xiàn)象，同時(shí)在圖像中灰度數(shù)量占比較少的微弱信息也會(huì)因?yàn)榛叶燃?jí)的拉伸而被合并[6]。圖像伽馬變換是另一種常用的空域圖像增強(qiáng)方法。伽馬校正可以有效地提升圖像的對(duì)比度，但是對(duì)于低照度圖像，普通的伽馬校正難以產(chǎn)生較好的效果。文獻(xiàn)[7]中提出了一種基于全局亮度平均值信息來(lái)對(duì)圖像局部區(qū)域進(jìn)行修正的伽馬校正圖像增強(qiáng)算法，該算法將圖像分為高通濾波和低通濾波兩個(gè)部分分別進(jìn)行操作，通過(guò)融合圖像的全局信息實(shí)現(xiàn)了圖像的增強(qiáng)。文獻(xiàn)[8]提出了一種改進(jìn)的伽馬校正函數(shù)，該方法通過(guò)圖像的局部信息自適應(yīng)地確定伽馬函數(shù)中的參數(shù)，有效地提高了圖像的對(duì)比度和亮度。

圖像的梯度反映了一副圖像的紋理信息，對(duì)圖像的梯度場(chǎng)進(jìn)行處理可以有效地提升圖像的質(zhì)量。目前針對(duì)于圖像梯度場(chǎng)進(jìn)行操作的圖像增強(qiáng)方法其步驟主要可分為以下三步：首先獲取原圖像的梯度場(chǎng)；其次對(duì)原梯度場(chǎng)進(jìn)行某種變換以獲得目標(biāo)梯度場(chǎng)；最后根據(jù)目標(biāo)梯度場(chǎng)進(jìn)行圖像重建獲得增強(qiáng)后的圖像[9]。在文獻(xiàn)[9]中，董麗麗、丁暢等通過(guò)對(duì)梯度場(chǎng)進(jìn)行局部均衡化，抑制了圖像高亮區(qū)域的擴(kuò)散，取得了不錯(cuò)的增強(qiáng)效果。文獻(xiàn)[10]中，趙文達(dá)等提出一種基于圖像梯度直方圖高斯規(guī)定化的紅外圖像增強(qiáng)算法，該算法通過(guò)對(duì)原圖像的梯度直方圖進(jìn)行高斯函數(shù)規(guī)定化，使紅外圖像的成像效果更加清晰。

本文在傳統(tǒng)伽馬校正的基礎(chǔ)上提出一種改進(jìn)的伽馬校正函數(shù)，該函數(shù)的參數(shù)融合了每個(gè)分塊圖像的全局信息。同時(shí)將改進(jìn)的伽馬變換作用于圖像的梯度域，獲取目標(biāo)梯度場(chǎng)后再進(jìn)行圖像的重建。在圖像重建階段，本文通過(guò)對(duì)圖像進(jìn)行分塊處理有效地減少了計(jì)算量，大大縮減了圖像的重建時(shí)間。另外，由于本文采用分塊重建策略重建了整個(gè)圖像，在算法具體實(shí)施時(shí)更具實(shí)用性。

2 算法原理

2.1 梯度場(chǎng)的閾值分割

本文在對(duì)圖像的梯度場(chǎng)進(jìn)行增強(qiáng)之前首先需要對(duì)其進(jìn)行梯度劃分，在圖像梯度場(chǎng)中，梯度值較小的像素點(diǎn)在所有的像素點(diǎn)中所占比例很高，體現(xiàn)在梯度直方圖中就是左高右低的類(lèi)似于“滑坡”的樣式，如圖1所示。針對(duì)于這種分布特點(diǎn)，文獻(xiàn)[11]提出一種新的梯度閾值劃分方法，即根據(jù)梯度直方圖的累計(jì)分布比例來(lái)劃分圖像梯度場(chǎng)的大梯度區(qū)間和小梯度區(qū)間。本文采用這種梯度劃分方法，在將圖像梯度場(chǎng)進(jìn)行劃分以后，對(duì)圖像的大梯度區(qū)間進(jìn)行適度增強(qiáng)操作；對(duì)于小梯度區(qū)間的梯度進(jìn)行大幅度的增強(qiáng)。以此來(lái)達(dá)到圖像梯度場(chǎng)共同增強(qiáng)的目的。

圖1 梯度直方圖

2.2 改進(jìn)的伽馬校正算法

伽馬校正早期被應(yīng)用于顯示器和掃描儀這一類(lèi)的成像設(shè)備上，應(yīng)用它的目的主要是對(duì)輸入圖像進(jìn)行校正。伽馬校正的具體形式見(jiàn)式(1)。

S(x,y)=C·R(x,y)γ

(1)

式中：R(x,y)表示輸入圖像的灰度值；S(x,y)表示經(jīng)過(guò)伽馬校正后的輸出圖像灰度值；C表示正值常數(shù)，為了將輸入灰度值和輸出灰度值歸一化為[0,1]的范圍，C值通常取1，因此式(1)可簡(jiǎn)化為S(x,y)=R(x,y)γ；式中的γ被稱(chēng)為伽馬校正系數(shù)，隨著其取值的不同圖像會(huì)產(chǎn)生不同的校正效果。

從圖2可以看出，當(dāng)γ>1時(shí)，函數(shù)輸出值會(huì)隨輸入值的增大而減小，具體表現(xiàn)為圖像的灰度值被壓縮，圖像變得更暗；反之，當(dāng)γ<1時(shí)，伽馬函數(shù)會(huì)增大原圖像的灰度值，圖像會(huì)變得更亮；當(dāng)γ=1時(shí)，伽馬函數(shù)將退化為恒等變換，不會(huì)對(duì)輸入圖像產(chǎn)生作用。因此為了達(dá)到增強(qiáng)梯度場(chǎng)的目的，本文所提出的改進(jìn)的伽馬函數(shù)主要在γ<1的系數(shù)范圍內(nèi)對(duì)原圖像梯度場(chǎng)進(jìn)行操作。

圖2 傳統(tǒng)伽馬校正函數(shù)圖像

傳統(tǒng)伽馬變換對(duì)于圖像的梯度場(chǎng)具有一定的增強(qiáng)作用，但其在加深梯度圖像背景的同時(shí)可能會(huì)拉低圖像中目標(biāo)物的梯度值[12]。此外，在傳統(tǒng)伽馬變換中，通常使用單個(gè)固定的伽馬校正系數(shù)對(duì)整張圖像進(jìn)行處理，這樣只能對(duì)目標(biāo)圖像進(jìn)行局部的增強(qiáng)。因此針對(duì)于以上兩點(diǎn)局限性，本文提出一種基于圖像梯度場(chǎng)局部信息的伽馬變換函數(shù)，見(jiàn)式(2)。

(2)

其中：

(3)

式中：Gout(x,y)表示伽馬函數(shù)輸出的圖像梯度值；Gin(x,y)表示輸入的圖像梯度值，在本文算法中，輸入梯度必須要?dú)w一化至[0,1]的區(qū)間范圍；Gmax表示原圖像的梯度最大值；T為2.1節(jié)所述的梯度圖像分割的梯度閾值，其數(shù)值表示的是梯度值而非梯度圖像分割的比例；θ和μ分別代表原圖像的梯度標(biāo)準(zhǔn)差和梯度平均值；β1和β2是調(diào)整系數(shù)。

在式(2)和式(3)中，β1和β2需要進(jìn)行設(shè)定，其余參數(shù)來(lái)源于梯度圖像的統(tǒng)計(jì)信息。圖3是基于圖像統(tǒng)計(jì)數(shù)據(jù)繪制的函數(shù)圖像。

圖3 改進(jìn)的伽馬函數(shù)圖像

從圖3可知，在不同的梯度分割區(qū)間內(nèi)，本文提出的伽馬函數(shù)對(duì)原圖像的梯度值采用了不同的放大策略。與大梯度值變換曲線相比，小梯度值變換曲線要更加陡峭一些。這體現(xiàn)出改進(jìn)的伽馬函數(shù)對(duì)大梯度區(qū)間進(jìn)行了一定程度的放大，而對(duì)于小梯度區(qū)間產(chǎn)生了更強(qiáng)的放大效果。這與2.1節(jié)中所提出的梯度分區(qū)間變換思路相符合。需要注意的是，由于大梯度區(qū)間的梯度值比較大，在經(jīng)過(guò)函數(shù)處理后會(huì)出現(xiàn)輸出值過(guò)大的情況，因此本文在算法實(shí)現(xiàn)階段對(duì)其輸出值進(jìn)行了限制，見(jiàn)式(4)。

Gout(x,y)=η·Gmax,Gout(x,y)>Gmax

(4)

對(duì)于低照度圖像，η取值一般在1～2.5范圍內(nèi)，本文取2。

2.3 圖像重建

對(duì)原圖像梯度場(chǎng)進(jìn)行增強(qiáng)處理后會(huì)得到目標(biāo)的梯度場(chǎng)G，因此需要尋找一個(gè)梯度場(chǎng)與目標(biāo)梯度場(chǎng) 最佳近似的灰度圖像U[13]。在數(shù)學(xué)上該過(guò)程等價(jià)于：在所有二維函數(shù)空間中，尋找一個(gè)梯度在最小二乘意義下與目標(biāo)梯度場(chǎng)G(x,y)最接近的函數(shù)U(x,y)，即求式(5)的泛函最小值[14]：

(5)

利用變分法[15]中的Euler-Lagrange方程對(duì)式(5)進(jìn)行化簡(jiǎn)整理最終可得式(6)：

(6)

可以看出式(6)就是泊松方程，因此對(duì)目標(biāo)梯度場(chǎng)的重建問(wèn)題可以轉(zhuǎn)化為求泊松方程的數(shù)值解問(wèn)題。而對(duì)泊松方程進(jìn)行數(shù)值求解實(shí)際上就是求解以下的線性方程組：

LU=divG

(7)

式中：L是Laplacian系數(shù)矩陣，該矩陣可以通過(guò)卷積方式構(gòu)建出來(lái)[16]；U是所求圖像矩陣的列向量形式。G是目標(biāo)梯度場(chǎng)的散度向量，求解式見(jiàn)式(8)。

divG=Gx(x,y)-G(x-1,y)+Gy(x,y)

-G(x,y-1)

(8)

上述方法雖然在算法的實(shí)現(xiàn)方面較為簡(jiǎn)單，但是仍然存在較大的問(wèn)題。以本文采用的512×512的圖像為例，在進(jìn)行圖像復(fù)原時(shí)，需要將目標(biāo)梯度場(chǎng)的散度矩陣按行拉伸成(262144,1)的列向量，因此與之相匹配的系數(shù)矩陣L的階數(shù)會(huì)迅速增加到(262144,262144)。過(guò)大的系數(shù)矩陣會(huì)造成數(shù)據(jù)存儲(chǔ)的困難，而且也很難對(duì)其進(jìn)行稀疏化。另外，在線性方程組的計(jì)算方面也會(huì)消耗大量的計(jì)算資源[10]。針對(duì)以上問(wèn)題，文獻(xiàn)[5]提出一種泊松方程快速解法，該方法是將二維的Laplacian系數(shù)矩陣分解為兩個(gè)正交方向上的1維矩陣，并對(duì)其進(jìn)行正弦變換，最后再利用Kronecker直積方法，將兩個(gè)1維向量組合為2維Laplacian系數(shù)矩陣。該方法解決了2維Laplacian系數(shù)矩陣的生成問(wèn)題，但系數(shù)矩陣的階數(shù)并未下降。在文獻(xiàn)[11]中，丁暢、董麗麗等在文獻(xiàn)[5]的基礎(chǔ)上提出一種矩陣變換法，該方法通過(guò)三角矩陣的相似變換對(duì)角化，大幅度地降低了系數(shù)矩陣的階數(shù)，計(jì)算速度顯著提高，但該方法在工程實(shí)現(xiàn)時(shí)需要多次矩陣變換稍顯復(fù)雜，因此希望能采用較為簡(jiǎn)單的策略解決這個(gè)問(wèn)題。

本文采用對(duì)圖像進(jìn)行分塊重建的方式降低計(jì)算量。以512×512的圖像為例，如果將原圖像劃分為4×4的圖像塊，則單塊圖像的系數(shù)矩陣會(huì)迅速下降為(16,16)，不但解決了系數(shù)矩陣存儲(chǔ)的問(wèn)題，還大大地縮減了圖像的重建時(shí)間。另外，這樣做避免了算法實(shí)現(xiàn)過(guò)程中的復(fù)雜變換，更具實(shí)際過(guò)程意義。

2.4 算法的總體流程

算法的總體流程見(jiàn)圖4。本文采取的方法在具體實(shí)現(xiàn)時(shí)存在以下幾點(diǎn)說(shuō)明：

圖4 算法的總體流程

(1)原圖像的梯度需要分成水平和豎直兩個(gè)方向分別求出。另外，因?yàn)閳D像的梯度值可能存在正負(fù)值及零值，因此需要用兩個(gè)矩陣分別保存兩幅梯度圖像的正負(fù)信息。

(2)在求取散度之前，需要將原梯度圖像的方向信息進(jìn)行還原。

(3)求取圖像的散度信息以后，需要將其按行壓縮為列向量，再通過(guò)2.3節(jié)方法重建后需要將圖像按照分塊的尺寸重構(gòu)成分塊圖像。

3 實(shí)驗(yàn)評(píng)價(jià)

3.1 本文方法效果的主觀評(píng)價(jià)

在圖像處理領(lǐng)域，對(duì)圖像具有增強(qiáng)作用的常用算法有直方圖均衡化和基于“視網(wǎng)膜大腦皮層理論”的Retinex算法(SSR、MSR、MSRCR)。因此本文在實(shí)驗(yàn)中選用文獻(xiàn)[5]算法和以上兩種算法與本文方法進(jìn)行比較。通過(guò)比較圖5至圖7發(fā)現(xiàn)，本文方法可以比較有效地提升MRI的對(duì)比度和細(xì)節(jié)輪廓。另外，直方圖均衡化等方法在處理圖像時(shí)對(duì)背景噪聲進(jìn)行了增強(qiáng)。相比較而言，本文方法噪聲幅值更小，在視覺(jué)體驗(yàn)上更加“干凈”。

3.2 本文方法效果的客觀評(píng)價(jià)

目前常用的圖像質(zhì)量評(píng)價(jià)方法有很多，本文采用信息熵[17]、圖像標(biāo)準(zhǔn)差和峰值信噪比作為算法的評(píng)價(jià)指標(biāo)。信息熵一般用來(lái)評(píng)價(jià)圖像中所含信息的豐富程度，表1是圖5至圖7中各類(lèi)算法的信息熵比較結(jié)果?？梢钥闯?，本文方法在一定程度上提升了原圖像的信息熵。由于直方圖均衡化等算法在圖像增強(qiáng)的過(guò)程中增大了圖像背景噪聲的強(qiáng)度，因此其較高的信息熵實(shí)際上體現(xiàn)了其過(guò)度增強(qiáng)的本質(zhì)。

表1 不同算法的信息熵對(duì)比

圖像的標(biāo)準(zhǔn)差是描述圖像細(xì)節(jié)的指標(biāo)，標(biāo)準(zhǔn)差值越大說(shuō)明圖像的細(xì)節(jié)表現(xiàn)得越清楚。表2的數(shù)據(jù)結(jié)果顯示了本文方法對(duì)于MRI具有較好的增強(qiáng)效果。

表2 不同算法的標(biāo)準(zhǔn)差對(duì)比

為了衡量在圖像增強(qiáng)過(guò)程中各類(lèi)算法是否同步地對(duì)噪聲進(jìn)行了增強(qiáng)，本文采用峰值信噪比(PSNR)來(lái)進(jìn)行各類(lèi)算法結(jié)果的評(píng)價(jià)。峰值信噪比在圖像處理領(lǐng)域中是一種非常常用的客觀評(píng)價(jià)指標(biāo)，它常用于評(píng)價(jià)圖像處理結(jié)果的失真程度。峰值信噪比越大，圖像處理結(jié)果失真越小，與原圖像的相似性越大。從表3中的數(shù)據(jù)可以看出，本文方法的處理結(jié)果與原圖像相似性更大，說(shuō)明本文方法只是對(duì)圖像特征進(jìn)行增強(qiáng)，而并未進(jìn)行圖像噪聲增強(qiáng)。

表3 不同算法的PSNR對(duì)比

3.3 本文算法速度提升評(píng)價(jià)

由2.3節(jié)可知，本文通過(guò)小分塊重建的方式對(duì)增強(qiáng)后的梯度場(chǎng)進(jìn)行圖像重建。在實(shí)驗(yàn)中發(fā)現(xiàn)，采用分塊的策略可以較大幅度地縮短算法的執(zhí)行時(shí)間，同時(shí)也能降低算法對(duì)內(nèi)存等計(jì)算資源的消耗。在圖5(e)、圖6(e)和圖7(e)中，本文采用4×4的分塊尺寸進(jìn)行重建，程序耗時(shí)分別是：5.24 s、5.38 s和5.49 s。由于本文采用的512×512分辨率，圖像尺寸過(guò)大，用傳統(tǒng)方法很難對(duì)其進(jìn)行整體重建，因此，另外選擇了一幅200×200分辨率的圖像來(lái)進(jìn)行對(duì)比測(cè)試。具體的計(jì)算時(shí)長(zhǎng)參照表4。

圖5 膝蓋部位的MRI增強(qiáng)效果圖及直方圖

圖6 腦部MRI增強(qiáng)效果圖及直方圖

從表4中可以看出，伴隨著分塊尺寸的縮小，本文算法的執(zhí)行速度得到較大的提升。這說(shuō)明本文采用的分塊重建的策略是成功的。在算法速度上基本滿(mǎn)足了實(shí)時(shí)性的要求。

表4 不同分塊尺寸的時(shí)長(zhǎng)對(duì)比(單位：s)

3.4 方法的不足之處

(1)本文方法在對(duì)梯度場(chǎng)進(jìn)行改進(jìn)的伽馬函數(shù)處理時(shí)未能完全實(shí)現(xiàn)參數(shù)自適應(yīng)，因此針對(duì)不同種類(lèi)的圖像需要調(diào)整不同的β1,β2參數(shù)值，但是過(guò)大的β1,β2可能會(huì)導(dǎo)致大、小梯度區(qū)間的變換曲線置于恒等變換曲線以下，即會(huì)對(duì)圖像梯度產(chǎn)生壓縮效果，因此針對(duì)于MRI圖像，如果參數(shù)取以下的范圍值可以獲得較好的效果：β1=6～9,β2=β1/2.5。

(2)通過(guò)觀察圖5至圖7發(fā)現(xiàn)，在采用本文方法進(jìn)行圖像分塊合并以后，圖像效果基本符合人眼的視覺(jué)要求。但如果在處理過(guò)程中分塊尺寸過(guò)大，經(jīng)過(guò)處理以后的圖像塊之間的灰度值會(huì)產(chǎn)生較大的差異，合并以后會(huì)產(chǎn)生比較明顯的塊狀效應(yīng)，這個(gè)弊端是本文采用的分塊策略導(dǎo)致的。

以上不足之處有待后續(xù)改進(jìn)。

4 結(jié)論

針對(duì)部分MRI圖像中存在的細(xì)節(jié)模糊不清、對(duì)比度低的問(wèn)題，本文提出一種基于改進(jìn)伽馬校正的圖像增強(qiáng)算法。該算法首先在原圖像梯度場(chǎng)上進(jìn)行改進(jìn)的伽馬校正函數(shù)處理，然后采用小分塊重建的方式對(duì)增強(qiáng)后的梯度場(chǎng)進(jìn)行重建，最后再合并小的圖像塊獲取增強(qiáng)后的圖像。通過(guò)與其他算法進(jìn)行比較可以看出，本文提出的方法對(duì)于MRI圖像具有較好的增強(qiáng)效果，且在運(yùn)算速度上也超過(guò)了傳統(tǒng)算法，對(duì)于本文采用的512×512分辨率圖片，在i3-10110U CPU和12 G內(nèi)存的硬件條件下，本文算法最快的處理時(shí)間約為5.24 s。此外，本算法也為傳統(tǒng)梯度圖像重建算法的改進(jìn)提供了一種新的思路。針對(duì)于本文方法的不足之處，有待后續(xù)進(jìn)行持續(xù)的改進(jìn)，以使其在具體應(yīng)用方面更加成熟。