基于物質(zhì)灰度正態(tài)分布的CT斷層圖像匹配插值

吳良武 侯建華

（大連測(cè)控技術(shù)研究所，大連 116013）

引言

斷層圖像插值是斷層圖像三維重建的關(guān)鍵技術(shù)之一，也是一個(gè)研究得比較多的問題。一般情況下，從CT/MRI等掃描設(shè)備獲得的斷層圖像，其層間距大于層內(nèi)像素間距。數(shù)據(jù)間的這種各向異性，將會(huì)影響醫(yī)學(xué)圖像分析及輔助治療系統(tǒng)中人體組織和器官的三維重建效果，也會(huì)給放射治療劑量計(jì)算帶來困難，因此有必要插值出各向同性的三維空間規(guī)則體數(shù)據(jù)［1］。

斷層醫(yī)學(xué)圖像插值方法可分為基于灰度的插值方法、基于形狀的插值方法和基于小波的插值方法?；谛螤畹牟逯到Y(jié)合從指定鄰域中提取出來的各種特征，用來計(jì)算插值點(diǎn)灰度。該方法通常用于三維重建面繪制技術(shù)，通過已知斷層圖像的形狀，直接構(gòu)造出中間插值圖像的輪廓，以方便顯示［2］。基于小波的插值算法，通過改變邊緣點(diǎn)對(duì)應(yīng)小波系數(shù)的位置和強(qiáng)度，對(duì)灰度和形狀同時(shí)插值，以充分模擬生理組織從一幅圖像變化到另一幅圖像的中間狀態(tài)［3－4］。

基于灰度的常用插值方法有線性插值、樣條插值和克立格（Kriging）插值等，從本質(zhì)上說都是計(jì)算采樣點(diǎn)灰度的加權(quán)平均，并且都根據(jù)插值點(diǎn)和采樣點(diǎn)之間的相對(duì)距離來調(diào)整權(quán)值的大小。其中，線性插值是樣條插值的特例。假定灰度在3個(gè)正交方向上都按樣條曲線的趨勢(shì)變化，用三元三次樣條張量積的形式來估計(jì)插值點(diǎn)的灰度。由于樣條插值預(yù)先假定了數(shù)據(jù)按某種樣條曲線變化，當(dāng)數(shù)據(jù)的實(shí)際變化規(guī)律與此有較大偏差時(shí)，則會(huì)產(chǎn)生不可忽略的錯(cuò)誤?？肆⒏癫逯祷跀?shù)據(jù)的空間分布規(guī)律，進(jìn)行統(tǒng)計(jì)意義上的最優(yōu)無偏估計(jì)?；叶炔逯捣椒ū举|(zhì)上是一種低通濾波，必然會(huì)對(duì)圖像中的高頻信息（如圖像的邊緣等）帶來畸變，因此這類處理方式產(chǎn)生的結(jié)果圖像普遍都存在邊緣模糊的問題。

針對(duì)傳統(tǒng)灰度插值方法所存在的問題，匹配插值的研究得到了人們的重視，比較典型的是Goshtasby等提出的基于對(duì)應(yīng)點(diǎn)匹配的插值方法［5－6］。該方法對(duì)于插值圖像中的每一個(gè)像素，根據(jù)插值點(diǎn)的位置選擇一系列的對(duì)應(yīng)點(diǎn)對(duì)，再通過計(jì)算對(duì)應(yīng)點(diǎn)對(duì)的特征相似性，選擇出最佳匹配點(diǎn)對(duì)，最后在最佳匹配點(diǎn)對(duì)間線性灰度插值。在Goshtasby方法的基礎(chǔ)上，有人提出了一種基于物質(zhì)分類的匹配插值方法［7］；首先計(jì)算各候選匹配點(diǎn)對(duì)為同種密度物質(zhì)的概率S，設(shè)置判定匹配點(diǎn)對(duì)為同種密度物質(zhì)的概率閾值，當(dāng)候選匹配點(diǎn)對(duì)的S值大于該閾值時(shí)，認(rèn)為點(diǎn)對(duì)中兩像素點(diǎn)分布于同種密度物質(zhì)，可作為匹配插值的進(jìn)一步候選點(diǎn)對(duì)，然后再使用Goshtasby方法進(jìn)行后續(xù)的匹配、插值等運(yùn)算。匹配插值將使得插值結(jié)果圖像的邊緣變得清晰。還有一種方法是：根據(jù)待插值圖像與其相鄰原始圖像的對(duì)應(yīng)像素的相關(guān)性，先對(duì)像素點(diǎn)進(jìn)行分類，然后再分別采用樣條插值、匹配點(diǎn)線性插值等不同方式，對(duì)不同分類點(diǎn)進(jìn)行插值，并對(duì)其進(jìn)行錯(cuò)誤校驗(yàn)［8］。

在當(dāng)前所看到的文獻(xiàn)中，當(dāng)在匹配點(diǎn)對(duì)間計(jì)算插值點(diǎn)的灰度值時(shí)，無一例外地都采用了以距離倒數(shù)作為權(quán)值的線性插值方法。筆者基于醫(yī)學(xué)CT斷層圖像中物質(zhì)灰度呈正態(tài)分布的特性，提出了一種新的插值方法：先分割CT斷層圖像，在分割過程中去除孤立噪聲點(diǎn)，然后篩選匹配點(diǎn)，使得候選匹配點(diǎn)對(duì)屬于同種物質(zhì)，為每個(gè)插值點(diǎn)搜索出各自的最佳匹配點(diǎn)對(duì)，最后在最佳匹配點(diǎn)對(duì)間，基于物質(zhì)灰度的正態(tài)分布曲線獲取插值權(quán)值，用來計(jì)算插值點(diǎn)灰度值。

1 CT斷層圖像分割

CT值表示對(duì)X光的吸收程度。當(dāng)一個(gè)像素處于兩種物質(zhì)的交界處，即存在部分體積效應(yīng)時(shí)，由于CT掃描設(shè)備分辨率有限，使得該像素的灰度值將位于兩種物質(zhì)的交叉處，如圖1所示。采用閾值分割法，CT切片灰度值由低到高，可分割為空氣、脂肪與軟組織、骨骼3種物質(zhì)，其中交叉處的像素屬于混合物質(zhì)。為解決部分體積效應(yīng)的分割問題，吳良武等提出了一種基于鄰域的骨骼確定分割算法［9］。

圖1 CT切片灰度分布Fig.1 Intensity distribution of CT slices

對(duì)于單個(gè)像素點(diǎn)，如果在切片內(nèi)，其8-鄰接像素都為其他物質(zhì)，而且在三維空間中，其6-鄰接像素也都為其他物質(zhì)，該像素將被認(rèn)為是孤立的噪聲點(diǎn)，被標(biāo)定為鄰域中權(quán)值最大的其他物質(zhì)。在切片內(nèi)，像素被看成小的正方形，8-鄰接表示像素間邊相連或者角相連。在三維空間中，像素被看成為體素，即小的立方體或長(zhǎng)方體，6-鄰接表示兩體素間面相連。

CT圖像分割后，根據(jù)文獻(xiàn)［5，7］的方法，設(shè)定匹配窗口寬度，在匹配窗口中選取屬于同種物質(zhì)的候選匹配點(diǎn)對(duì)，分別以匹配點(diǎn)間灰度值的相似度、兩點(diǎn)為中心區(qū)域的灰度變化的相似度及兩點(diǎn)水平坐標(biāo)接近程度等信息為權(quán)值，從候選匹配點(diǎn)對(duì)中為每個(gè)插值點(diǎn)搜索出各自的最佳匹配點(diǎn)對(duì)。

2 基于物質(zhì)灰度呈正態(tài)分布特性的插值方法

把連續(xù)CT斷層切片的灰度值看作隨機(jī)量，直方圖可以看作其概率密度函數(shù)（PDF）的估計(jì)p（z）［10］。取一組連續(xù)50層的人體頭部 CT切片，圖2是該組切片的脂肪和軟組織灰度直方圖、骨骼灰度直方圖，可以看出直方圖的形狀都相似于正態(tài)分布曲線。筆者對(duì)另外3組醫(yī)學(xué)CT斷層切片也分別繪制了脂肪和軟組織灰度直方圖、骨骼灰度直方圖，發(fā)現(xiàn)它們的形狀也都與正態(tài)分布曲線相似。

圖2 灰度直方圖。（a）脂肪和軟組織；（b）骨骼Fig.2 Intensity histogram.（a）Fat and soft tissues；（b）Bones

式中，μ和σ2分別表示均值和方差。

一般來說，如果一個(gè)像素全部由同種物質(zhì)組成，則它的CT值應(yīng)接近于某 CT0值。交叉像素由于部分體積效應(yīng)，其CT值將很難穩(wěn)定在 CT0值附近，使得CT斷層切片中同種物質(zhì)的灰度分布可能呈現(xiàn)出正態(tài)分布的特性。人體頭部存在大塊的脂肪和軟組織，其交叉像素所占比例小，方差σ2值小，分布圖看起來尖銳；人體頭部骨骼由多個(gè)小塊組成，交叉像素所占比例大，方差σ2值大，分布圖看起來平坦。圖2所示的灰度直方圖正符合這樣的解釋。

基于上述分析，可認(rèn)為在醫(yī)學(xué)CT斷層切片中，物質(zhì)灰度呈正態(tài)分布特性。

灰度直方圖擬合為正態(tài)分布曲線，能夠去除細(xì)節(jié)，給出其光順分布趨勢(shì)。例如，擬合后，圖2（b）中明顯的曲線細(xì)節(jié)將會(huì)被光順。圖3顯示了另一組脂肪和軟組織的灰度直方圖，CT掃描設(shè)備間隔取值，擬合曲線能去除取值間隔，給出光順分布曲線。使用共軛梯度爬山法，將能使擬合的均分誤差降至最?。?0］。

引入約束條件“匹配點(diǎn)對(duì)間灰度值的概率分布與該物質(zhì)的灰度正態(tài)分布一致”，由該物質(zhì)的灰度正態(tài)分布定義為［11］正態(tài)分布曲線來獲取計(jì)算插值點(diǎn)的權(quán)值，利用這種方法計(jì)算出的插值點(diǎn)灰度值將符合該物質(zhì)的正態(tài)分布特性。圖4是根據(jù)圖2（b）繪制的骨骼灰度積分分布曲線，由此容易得到用于計(jì)算新插值點(diǎn)灰度值的權(quán)值［11］。

圖3 另一組脂肪和軟組織的灰度直方圖Fig.3 Another intensity histogram of fat and soft tissues

圖4 骨骼灰度積分分布Fig.4 Intensity in tegraldistribution of bones

傳統(tǒng)的同種物質(zhì)匹配點(diǎn)對(duì)間以距離倒數(shù)作為權(quán)值的線性插值［5，7］，實(shí)際上隱含了這樣的假設(shè)條件：在匹配點(diǎn)對(duì)的連線上，灰度值呈線性規(guī)律變化。線性變化的概率密度在統(tǒng)計(jì)學(xué)中是均勻分布的，其積分分布曲線為直線［11］。均勻分布與物質(zhì)灰度呈正態(tài)分布特性是相矛盾的。如果匹配點(diǎn)對(duì)的灰度值相差不大，積分分布曲線在小區(qū)段間將可近似為直線。例如，它們都位于正態(tài)分布的頂部，或者都位于遠(yuǎn)離正態(tài)分布頂部的同一側(cè)，那么這兩點(diǎn)間可近似為均勻分布，此時(shí)使用線性插值與使用正態(tài)分布插值所得的結(jié)果基本沒有區(qū)別。不失一般性，假設(shè)匹配點(diǎn)對(duì)的灰度值相差較大，一個(gè)位于正態(tài)分布頂部，另一個(gè)位于遠(yuǎn)離頂部的一側(cè)（如圖5所示），此時(shí)計(jì)算出來的結(jié)果u和v將存在較大的差異。

3 測(cè)試

圖5 插值點(diǎn)值計(jì)算Fig.5 Calculating the value of the interpolating pixel

以上節(jié)提及的那組50層連續(xù)人體頭部CT切片作為測(cè)試用例。該組50層切片是從頭頂部位向下順序生成，切片間距為1.84 mm，圖6列出了其中的前10幅連續(xù)切片，圖中采用字母表示切片序號(hào)，測(cè)試結(jié)論將主要針對(duì)這10幅連續(xù)切片給出。去除切片中與人體頭部無關(guān)的背景固定物，定義灰度值大于135的像素為確定骨骼，灰度值介于65～115之間的像素為確定脂肪與軟組織，灰度值小于45的像素為確定空氣，灰度值介于115～135之間及介于45～65之間的像素為不確定像素。

每次測(cè)試選取3幅連續(xù)切片，從前、后切片中插值新切片，然后與中間原始切片做比對(duì)，計(jì)算插值與真實(shí)值間的偏差。

圖6 連續(xù)CT頭部切片。（a）～（j）為從頭頂部位開始的切片序列Fig.6 Continuous CT slices of the head.（a）～（j）The sequences of slices beginning at the top of the head

CT切片分割后，按照文獻(xiàn)［5，7］的方法，先為每一插值點(diǎn)搜索屬于同種物質(zhì)的最優(yōu)匹配點(diǎn)對(duì)，然后在最優(yōu)匹配點(diǎn)對(duì)間，分別使用線性插值與正態(tài)分布插值，計(jì)算插值點(diǎn)的灰度值。只討論兩種方法插值結(jié)果灰度值差大于2的插值點(diǎn)，對(duì)它們作標(biāo)記，采用標(biāo)記點(diǎn)的均偏差對(duì)插值結(jié)果的優(yōu)劣進(jìn)行評(píng)價(jià)。以N表示每次測(cè)試的標(biāo)記點(diǎn)總數(shù)，v表示計(jì)算出的標(biāo)記點(diǎn)插值，v0表示標(biāo)記點(diǎn)真實(shí)值，每次測(cè)試的均偏差e定義為

測(cè)試分兩組進(jìn)行。第1組測(cè)試在確定像素間進(jìn)行，也就是說，候選匹配點(diǎn)對(duì)都為同種物質(zhì)的確定像素。第2組測(cè)試規(guī)定只要有一個(gè)像素為確定像素，另一個(gè)像素如果為該物質(zhì)的不確定像素，則它們也可以作為候選匹配點(diǎn)對(duì)。這樣，第2組測(cè)試中將包含較多的匹配點(diǎn)對(duì)，一個(gè)為確定像素，另一個(gè)為不確定像素。

4 結(jié)果

測(cè)試結(jié)果如表1和圖7所示。表1列出了中間比對(duì)切片序號(hào)為（b）～（h）的7次測(cè)量數(shù)值結(jié)果。EN1列列出了第1組測(cè)試使用正態(tài)分布插值的標(biāo)記點(diǎn)均偏差值，EL1列列出了第1組測(cè)試使用線性插值的標(biāo)記點(diǎn)均偏差值；EN2列列出了第2組測(cè)試使用正態(tài)分布插值的標(biāo)記點(diǎn)均偏差值，EL2列列出了第2組測(cè)試使用線性插值的標(biāo)記點(diǎn)均偏差值。

將表1的測(cè)試結(jié)果繪制于圖7中。圖7的橫坐標(biāo)表示每次測(cè)試的中間比對(duì)切片序號(hào)，其中的橫坐標(biāo)2對(duì)應(yīng)切片（b），橫坐標(biāo)3對(duì)應(yīng)切片（c），依此類推，縱坐標(biāo)表示每次測(cè)試標(biāo)記點(diǎn)的均偏差 e。50幅連續(xù)切片共進(jìn)行48次測(cè)試，圖7列出了其中的前20次測(cè)試結(jié)果，繪制出了第1組線性插值和正態(tài)分布插值的均偏差曲線，以及第2組線性插值和正態(tài)分布插值的均偏差曲線。

表1 測(cè)試結(jié)果Tab.1 Experimental results

圖7 均偏差比較Fig.7 Comparison of average errors

先比較線性插值與正態(tài)分布插值的測(cè)試結(jié)果。從圖7可以看出，在前4次測(cè)試中，用正態(tài)分布插值方法得到的兩組均偏差明顯都比用線性插值方法的小。表1的數(shù)據(jù)也說明了這一點(diǎn)：在前4次測(cè)試中，第1組正態(tài)分布插值得到的均偏差比線性插值得到的均偏差平均減小了21.5％，第2組正態(tài)分布得到的均偏差比線性插值得到的均偏差平均更是減小了29.55％。不過，后續(xù)的測(cè)試結(jié)果則很難看出其中的明顯差別。圖7未列出的后28次測(cè)試結(jié)果也很難看出明顯差別。

查看圖6，從（a）～（f）的切片位于頭部的頂部，用于測(cè)試的前、后兩幅切片對(duì)應(yīng)像素的灰度值差別比較大，使得前4次測(cè)試所搜索到的最佳匹配點(diǎn)對(duì)中，存在灰度值差大的匹配點(diǎn)對(duì)較多。例如，在使用切片（a）和切片（c）插值的測(cè)試中，可以觀察到切片（c）中很多骨骼像素在切片（a）中的對(duì)應(yīng)匹配點(diǎn)的灰度值都比較小，導(dǎo)致最佳匹配點(diǎn)間差值大的插值點(diǎn)增多，這種現(xiàn)象直到中間切片序號(hào)為（e）的測(cè)試中都有存在。但在后續(xù)的測(cè)試中，前、后切片相似度比較高，如切片（f）和切片（h）之間的相似性要高很多，導(dǎo)致搜索到的最佳匹配點(diǎn)間差值并沒有那么明顯。在前4次測(cè)試中，因?yàn)檎龖B(tài)分布插值方法改進(jìn)了線性插值方法處理匹配點(diǎn)間差值大的情況，使得標(biāo)記點(diǎn)的均偏差變小，后續(xù)測(cè)試均偏差沒有明顯變化，是因?yàn)楹竺娴那啊⒑笄衅嗨贫缺容^高，導(dǎo)致搜索到的最佳匹配點(diǎn)間差值大的插值點(diǎn)少。此時(shí)，如上節(jié)所述，兩種插值方法得到的結(jié)果沒有大的差別。

再比較第1組與第2組的測(cè)試結(jié)果。從圖7可以看出，在前4次測(cè)試中，線性插值在第2組測(cè)試中得到的標(biāo)記點(diǎn)均偏差明顯比第1組測(cè)試增大，從表1可計(jì)算出平均增大了14.43％，正態(tài)分布插值方法的標(biāo)記點(diǎn)均偏差在兩組測(cè)試中則基本保持不變。這說明，隨著匹配點(diǎn)間差值大的插值點(diǎn)增多，線性插值方法的標(biāo)記點(diǎn)均偏差也在變大，而正態(tài)分布插值方法的標(biāo)記點(diǎn)均偏差基本保持不變，表明正態(tài)分布插值方法與實(shí)際情況更相符。

上述兩種測(cè)試分析結(jié)果都說明，基于正態(tài)分布的插值方法改進(jìn)了當(dāng)前最佳匹配點(diǎn)對(duì)間普遍采用的線性插值方法，正確修正了其中的匹配點(diǎn)間差值大的插值點(diǎn)的插值，從而整體上提高了插值效果。

5 結(jié)論

本研究基于CT斷層切片同種物質(zhì)灰度呈正態(tài)分布的特性，引入了約束條件“匹配點(diǎn)對(duì)間灰度值的概率分布與該物質(zhì)的灰度正態(tài)分布一致”，在最佳匹配點(diǎn)對(duì)間，基于物質(zhì)的灰度正態(tài)分布獲取權(quán)值，計(jì)算插值點(diǎn)的灰度值。測(cè)試結(jié)果表明，這種方法改進(jìn)了當(dāng)前匹配點(diǎn)對(duì)間普遍采用的線性插值，正確修正了其中的匹配點(diǎn)間差值大的插值點(diǎn)的插值，從而整體上提高了插值效果。對(duì)于其他類型的醫(yī)學(xué)斷層切片匹配插值，該插值方法應(yīng)該也具有一定的參考應(yīng)用價(jià)值。

［1］陳武凡，秦安，江少峰，等.醫(yī)學(xué)圖像分析的現(xiàn)狀與展望［J］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，2008，27（2）：175－181.

［2］GreveraGJ， Udupa JK. Shape － based interpolation of multidimensional gray-level images［J］.IEEE Trans on Medical Imaging，1996，15（6）：881－892.

［3］黃海赟，戚飛虎，陳劍，等.基于小波的醫(yī)學(xué)圖像插值［J］.自動(dòng)化學(xué)報(bào)，2002，28（5）：722－728.

［4］胡俊峰，唐鶴云，錢建生.基于小波變換醫(yī)學(xué)圖像融合算法的對(duì)比分析［J］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，2011，30（2）：196－205.

［5］Goshtasby A，Turner DA.Matching of tomographic slices for interpolation［J］. IEEE Transactions on Medical Imaging，1992，11（4）：507－516.

［6］繆曉和，鄧元木，黃斐增，等.基于對(duì)應(yīng)點(diǎn)匹配的斷層圖象三維插值方法［J］.中國醫(yī)學(xué)物理學(xué)雜志，2000，17（1）：14－16.

［7］張勇，歐宗瑛，秦緒佳，等.基于物質(zhì)分類的三維空間斷層圖像匹配插值［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2002，14（7）：559－663.

［8］田沄，齊敏，衛(wèi)旭芳，等.基于像素分類的醫(yī)學(xué)圖像層間插值［J］.中國圖象圖形學(xué)報(bào)，2008，13（9）：1655－1660.

［9］吳良武，侯建華，張勇，等.用鄰域運(yùn)算從 CT圖像中分割骨骼［J］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，2003，22（3）：199－202.

［10］岡薩雷斯，著.阮秋琦，譯.數(shù)字圖像處理（第2版）［M］.北京，電子工業(yè)出版社，2005，489－493.

［11］馬洪寬，張華隆，主編.概率論與數(shù)理統(tǒng)計(jì)［M］.上海：上海交通大學(xué)出版社，2005.