王德兵

(安徽工貿(mào)職業(yè)技術(shù)學(xué)院計(jì)算機(jī)信息工程系，安徽 淮南 232007)

醫(yī)學(xué)圖像配準(zhǔn)一直是研究者普遍關(guān)注的問題?，F(xiàn)在，已開發(fā)出多種方法可用于醫(yī)學(xué)圖像配準(zhǔn)。其中，基于標(biāo)志點(diǎn)的圖像配準(zhǔn)在現(xiàn)實(shí)中最為常用，如在諸多商業(yè)圖像導(dǎo)航系統(tǒng)中廣泛使用此技術(shù)，而基于標(biāo)志點(diǎn)的圖像配準(zhǔn)涉及不同圖像(或物理空間)中對應(yīng)點(diǎn)坐標(biāo)的確定，以及使用這些對應(yīng)點(diǎn)開展的幾何變換估計(jì)，這些標(biāo)記點(diǎn)可以是內(nèi)在標(biāo)志或外部標(biāo)志。內(nèi)在標(biāo)志一般來源于自然存在的特征，如解剖標(biāo)志點(diǎn)，而外部標(biāo)志則為人工選擇標(biāo)志。基準(zhǔn)點(diǎn)是基于標(biāo)記的顱腦圖像配準(zhǔn)中最重要的參考指標(biāo)之一，而顱腦圖像中基準(zhǔn)點(diǎn)位置的估算對于成功開展基于標(biāo)記的圖像配準(zhǔn)十分重要。現(xiàn)有的商業(yè)計(jì)算機(jī)輔助手術(shù)軟件系統(tǒng)要么不具備自動定位基準(zhǔn)點(diǎn)的功能，要么聲稱具有此功能，但實(shí)際上并非如此。比如，Stealth Station(手術(shù)導(dǎo)航系統(tǒng))就不具備此功能，史塞克(Stryker)在系統(tǒng)界面上設(shè)置“自動發(fā)現(xiàn)(Find automatically)”按鍵并聲稱具備此功能，然而實(shí)際顯示所有的自動發(fā)現(xiàn)功能均不能正常工作，史塞克醫(yī)療器械有限公司承認(rèn)確實(shí)不能使用這些功能，這就意味著基準(zhǔn)標(biāo)記仍需用戶手動選擇，用戶必須觀察所有圖像切面并使用鼠標(biāo)選擇基準(zhǔn)中心用于配準(zhǔn)。在手動選擇配準(zhǔn)過程中，圖像強(qiáng)度、對比度、切面間距，甚至是圖像中不完整基準(zhǔn)點(diǎn)均會使不同用戶生成不同的鼠標(biāo)定位結(jié)果。因此人工選擇具有較大缺陷，并可能造成誤差。自動基準(zhǔn)點(diǎn)檢測是圖像配準(zhǔn)的首要步驟，近年來研究者關(guān)于這一主題已開展了一些研究。文獻(xiàn)[1]提出了一種基于強(qiáng)度的搜索和分類的自動標(biāo)記定位技術(shù)，然而，這一方法需要關(guān)于標(biāo)記大小和形狀的先驗(yàn)信息，且僅適用于圓柱形標(biāo)記。文獻(xiàn)[2]開發(fā)了基于超聲波模塊的方法用于檢測并定位皮下基準(zhǔn)點(diǎn)，但該方法需使用附加硬件設(shè)備才能完成基準(zhǔn)點(diǎn)的定準(zhǔn)。文獻(xiàn)[3]也提出了相應(yīng)的方法，但要求基準(zhǔn)點(diǎn)僅限于最佳尺寸大小。文獻(xiàn)[4]提出了種子和條件膨脹方法，而該方法則需要高清晰圖像。其他技術(shù)需要則需要增加硬件設(shè)備，如關(guān)節(jié)臂[5-6]、光學(xué)三角系統(tǒng)[7]、磁場數(shù)字轉(zhuǎn)換器[8]等。本文充分利用圖像空間的強(qiáng)度信息，提出了一種基于模板的基準(zhǔn)點(diǎn)檢測算法和決策系統(tǒng)，完成基準(zhǔn)點(diǎn)的定準(zhǔn)，克服了上述方法的缺陷。

1 檢測系統(tǒng)結(jié)構(gòu)

基準(zhǔn)點(diǎn)檢測過程如圖1所示。

圖1 基準(zhǔn)點(diǎn)檢測過程

通常情況下，使用掃描儀器獲取的圖像為軸位圖像切面。切面間距(即兩個連續(xù)切面間的距離)可根據(jù)模態(tài)和醫(yī)生的要求固定，最常見的切面間距為2mm，一個圖像序列可含有任意數(shù)目的圖像切面。在大多數(shù)情況下，全腦圖像序列含有75張切面，部分大腦圖像序列含有25或52張切面，集中于腫瘤區(qū)域。一個圖像序列可創(chuàng)建三維圖像，即三維矩陣M(x，y，z)，M(x，y，z)的每個元素均是切面z中像素(x，y)的強(qiáng)度值。在大多數(shù)情況下，0≤x，y≤255， 0≤z≤25(或者0≤z≤50，0≤z≤75)。從M(x，y，z)中，提取了三種不同的圖像視圖：即軸位、冠狀位和矢狀位視圖。軸位視圖實(shí)際是輸入圖像序列，并對三個視圖分別進(jìn)行處理。在這一階段使用圖像處理技術(shù)是為了獲取皮膚和基準(zhǔn)點(diǎn)清晰的輪廓，如果一段輪廓的曲率具有基準(zhǔn)點(diǎn)的特征，這一段輪廓被選為候選基準(zhǔn)點(diǎn)。從三個視圖中獲取三個候選集，作為基于模板的決策系統(tǒng)的輸入，再使用該系統(tǒng)測試軸位視圖中各個候選基準(zhǔn)點(diǎn)，以確定在相同的三維位置是否也有候選基準(zhǔn)點(diǎn)。如果有，認(rèn)為這個候選點(diǎn)是基準(zhǔn)點(diǎn)，否則，這個候選點(diǎn)是邊緣檢測的誤差。誤差可能有諸多原因?qū)е?，例如，圖像中不完整的基準(zhǔn)點(diǎn)，或一段輪廓是從某個角度看起來像是基準(zhǔn)點(diǎn)但卻是皮膚的一部分。

2 方法

本研究關(guān)鍵是提出了從三維空間檢測基于模板的自動檢測基準(zhǔn)點(diǎn)的算法，作為決策系統(tǒng)的核心，解決了基準(zhǔn)點(diǎn)檢測自動定準(zhǔn)問題。通過大腦邊緣圖的構(gòu)建，獲取大腦輪廓；根據(jù)大腦輪廓中角點(diǎn)的不同曲率確定候選基準(zhǔn)點(diǎn)；根據(jù)上述方法獲取軸位、冠狀位和矢狀位三個視圖的三個侯選集，作為決策系統(tǒng)的輸入。本部分分別詳細(xì)介紹了算法中使用到的模板，邊緣檢測方法，基于曲率的候選基準(zhǔn)點(diǎn)檢測方法和決策系統(tǒng)中基于模板的三維空間檢測算法。

2.1 基準(zhǔn)點(diǎn)模板

圖2為研究中使用的基準(zhǔn)標(biāo)記模型。如將該模型放在頭骨上進(jìn)行掃描，可得到圖3中所示的9個剖視圖。基于剖視圖，掃描平面會切割基準(zhǔn)標(biāo)記。例如，如果剖切面與z軸垂直，可得到剖視圖3(g)。如果剖切面與x或y軸垂直，可分別得到剖視圖3(h)和3(i)，其他6個剖視圖是在剖切面與任一軸均不垂直的情況下獲得的。

圖2 基準(zhǔn)標(biāo)記模型

圖3 剖切面在不同位置時(shí)的剖視圖

將上述9個視圖分為四類：視圖(e)、 (f) 和 (g)屬于第0類，視圖(h)屬于第1類，視圖(i)屬于第2類，所余視圖屬于第3類。本文方法僅使用第0、1、 2類視圖。

對于圖4中的三維圖像，當(dāng)對含有完整基準(zhǔn)標(biāo)記的連續(xù)圖像序列開展邊緣檢測時(shí)，有以下兩種可能(分別使用以下模板進(jìn)行說明)。

圖4 大腦頂部基準(zhǔn)標(biāo)記的多平面視圖

模板1：至少有第0類的剖視圖，如圖3中的(a)和(b)。圖4中第一行為基準(zhǔn)標(biāo)記的掃描圖，第二行為所檢測到的這些標(biāo)記的邊緣。所檢測到的邊緣為兩個不被圖像其他邊緣完全包圍的同心圓，這就保證了圖像只能是頭骨上設(shè)置的基準(zhǔn)標(biāo)記。

圖5 基準(zhǔn)標(biāo)記軸位圖及三幅連續(xù)切面中檢測到的邊緣

模板2：依次至少有一幅第0類剖視圖，至少有一幅第2類剖視圖和至少有一幅第1類剖視圖組成的圖像序列。在此圖像序列前后，均無基準(zhǔn)標(biāo)記。在第一類剖視圖中，邊緣含有單峰，稱之為“邊緣類型1”；在第2類剖視圖中，邊緣含有雙峰，稱為“邊緣類型2”；另外還有一種情況，邊緣中無峰，稱之為“邊緣類型0”。這樣，根據(jù)基準(zhǔn)標(biāo)記圖像序列中邊緣上的峰的數(shù)量，可分別從圖6和圖7中獲得字符串“0111112221000” 和“01122222211100”。對于連續(xù)的相同數(shù)字，僅保留一個，將多余的都刪除，以此方式可將圖6和圖7的字符串均收縮為“01210”。

圖6 圖3中基準(zhǔn)標(biāo)記的矢狀位視圖，13幅連續(xù)切面中檢測到的邊緣，以及邊緣上峰的數(shù)目

圖7 圖3中基準(zhǔn)標(biāo)記的冠狀位視圖，14幅連續(xù)切面中檢測到的邊緣，以及邊緣上峰的數(shù)目

以上模板1和2中，同心圓和峰值均在僅有一個基準(zhǔn)點(diǎn)的區(qū)域內(nèi)進(jìn)行檢測，研究中選擇基準(zhǔn)點(diǎn)最大可能尺寸為12mm。

2.2 邊緣圖構(gòu)建

使用基于比值的閾值法[9]在原始圖像中將前景對象從背景中分割出來，對于最佳閾值，兩個區(qū)域像素?cái)?shù)的比值會發(fā)生明顯變化。從所定位的潛在閾值集中，可自動確定最終閾值的子集，第一個閾值將用于前景/背景分割。使用比值曲線進(jìn)行分割的步驟如下：

生成直方圖并使用高斯濾波器對其進(jìn)行平滑處理，以抑制較小的強(qiáng)度變化。

計(jì)算比值曲線的二次倒數(shù)，獲得歸一變化率曲線，然后將曲線平滑處理，以消除較小的峰值和谷值，所得到的最大和最小值的位置決定了候選閾值的值，示例圖像的分割結(jié)果如圖8所示。

(a) 原始圖

(b) 分割圖圖8 前景/背景分割

在分割后的前景中，使用遞歸形態(tài)侵蝕、膨脹和灰度重建技術(shù)融合較小間隙和消除多余噪聲(常見于病例圖像中)，然后使用Canny邊緣檢測技術(shù)以獲取邊緣[10](圖9a)。

在檢測到的所有邊緣中，使用射線搜索程序選擇大腦的輪廓：1)對每個邊緣目標(biāo)的像素進(jìn)行編號并選擇上半部分作為候選2)從圖像四個角向中心設(shè)置四條射線3)與任一射線相交的第一條邊緣即大腦的輪廓。在開展此程序前，需使用高斯平滑濾波器消除較小的物體和噪聲(圖9b)。

(a) 初始邊緣圖

(b) 輪廓選擇圖9 輪廓選擇

2.3 基于曲率的候選基準(zhǔn)點(diǎn)檢測

大腦輪廓識別后，檢測輪廓內(nèi)的角點(diǎn)。單尺度特征檢測不能同時(shí)檢測細(xì)微和粗略特征，而多尺度特征檢測[11]可解決這一問題。研究中使用自適應(yīng)曲率閾值，而非單一的全局閾值。文中，曲率定義為：

在固定低尺度下計(jì)算曲率，以保留所有真實(shí)角點(diǎn)。曲率的所有局部最大值都被視為候選角點(diǎn)，包括錯誤角點(diǎn)，使用自適應(yīng)局部閾值和角點(diǎn)角度移除初始候選角點(diǎn)中的錯誤角點(diǎn)和邊界噪聲。角點(diǎn)檢測結(jié)果如圖10所示。

圖10 角點(diǎn)檢測

根據(jù)角點(diǎn)位置，使用K-均值聚類將所選擇的角點(diǎn)聚類至若干個類別中，再進(jìn)行多項(xiàng)式擬合，搜索可將數(shù)據(jù)擬合至各個類別的三級多項(xiàng)式P(x)的系數(shù)，繪制各組的擬合曲線，如圖11所示。選擇各條曲線的中心點(diǎn)，投射至大腦輪廓上，所得到的點(diǎn)即為基準(zhǔn)標(biāo)記的中心(見圖12)。

圖11 曲線擬合

圖12 最終結(jié)果

2.4 基于模板的三維空間基準(zhǔn)點(diǎn)檢測

從三個視圖獲取三個候選基準(zhǔn)點(diǎn)集合后，使用決策系統(tǒng)確定基準(zhǔn)點(diǎn)的位置。本節(jié)給出了決策系統(tǒng)及基準(zhǔn)點(diǎn)檢測算法的規(guī)則。對于從DICOM圖像序列重建的三維圖像“IMG”(包含N切面)，令I(lǐng)(a，i)，I(s，j)，和I(c，k)表示軸位切面圖像i，矢狀位切面圖像j，和冠狀位切面圖像k。

三維空間基準(zhǔn)點(diǎn)檢測算法描述如下：

1. 定義Ma，Ms和Mc為與“IMG”具有相同尺寸的矩陣，且所有元素值設(shè)為0。

2. 令i=1。

3. 對I(a，i)開展邊緣檢測。

4.開展基于曲率的候選基準(zhǔn)點(diǎn)檢測。如果“模板1”存在，令Ma(x，y)=1，其中 (x，y)是同心圓中心坐標(biāo)。

5.遍歷整個大腦邊界。對于所遇到的每個“邊緣類型2”，對I(a，i-1)，I(a，i-2)，…， 和I(a，i+1)，I(a，i+2)等開展邊緣檢測，在此圖像序列中相同大腦邊界位置尋找峰數(shù)字符串，如果有“模板2”，令Ma(x，y)=2，其中 (x，y)是兩個峰之間谷值坐標(biāo)。

6.如果既沒有“模板1”也沒有“模板2”，并且不存在i=N，則令i=i+1，返回第一步。

7. 對Ms中所有值為1的元素，即Ms(x，y，z)=1，進(jìn)行以下操作：

a) 令j=x，使用I(s，j)替換I(a，i)且使用Ms替換Ma，重復(fù)第3～5步操作；

b)令k=y，使用I(c，j)替換I(a，i)且使用Mc替換Ma，重復(fù)第3～5步操作。

8. 如操作7(a) 和 7(b)中均有“模板2”，則在坐標(biāo)(x，y，z)處有基準(zhǔn)標(biāo)記。

9. 對Ms中所有值為2的元素，即Ms(x，y，z)，=2，開展以下操作：

a) 令j=x，使用I(s，j)替換I(a，i)且使用Ms替換Ma，重復(fù)第3～5步操作；

b)令k=y，使用I(c，j)替換I(a，i)且使用Mc替換Ma，重復(fù)第3～5步操作。

10.如操作9(a) 和 9(b)中有“模板1”或“模板2”，則在坐標(biāo)(x，y，z)處有基準(zhǔn)標(biāo)記。

發(fā)現(xiàn)檢測邊緣以找到“模板2”的峰值并非易事，通過多種邊緣檢測濾波器如“Sobel”，“Laplacian”, “Prewiit”, “Zero-cross”, “Canny”等檢測發(fā)現(xiàn)，其中使用“Laplacian”濾波器獲取結(jié)果最優(yōu)，所有結(jié)果均是基于此邊緣檢測算法。同時(shí)，盡管“Laplacian”在邊緣檢測中仍然不夠完善，但通過考慮連續(xù)圖像序列而非單幅圖像，仍可以成功發(fā)現(xiàn)“模板2”并最終找到基準(zhǔn)標(biāo)記。

3 實(shí)驗(yàn)結(jié)果

盡管模板匹配方法在發(fā)現(xiàn)和定位標(biāo)記中非常有用，但使用該方法時(shí)需將標(biāo)記的基準(zhǔn)點(diǎn)定義為其形心。該方法首先將標(biāo)記從背景中分割出來，然后計(jì)算其強(qiáng)度加權(quán)形心。由于所使用的標(biāo)記在CT和MR圖像中較亮，在空氣背景上無法成像。因此，分割時(shí)將會遇到以下問題：使用閾值法將圖像分割為標(biāo)記和背景體元，和與區(qū)域生長相關(guān)的標(biāo)記體元。

因此本研究使用自適應(yīng)閾值法RaBAM[12-13]。Sahoo 和 Glasbey對現(xiàn)有的閾值算法進(jìn)行評閱和對比發(fā)現(xiàn)這些技術(shù)通常假設(shè)前景和背景體元來自兩個主要模式。因?yàn)闃?biāo)記的尺寸與圖像切面的厚度在同一數(shù)量級，許多含有標(biāo)記的體元會遭受部分容積效應(yīng)，即這些體元包含標(biāo)記和空氣的混合物。因此，圖像在標(biāo)記周圍區(qū)域的強(qiáng)度直方圖并沒有兩個主要模式。通常，強(qiáng)度直方圖有一個對應(yīng)背景噪聲的主要模式和對應(yīng)前景標(biāo)記的廣泛分布區(qū)域。如果含有組織或者CT圖中的標(biāo)志桿，背景模式也有可能是分散的。對于算法的第一部分，我們發(fā)現(xiàn)RaBAM算法可以勝任。

本研究使用六位病人的42套醫(yī)學(xué)圖像(每人一套CT和六套MR圖像，正向和反向讀出梯度T1，PD，T2)開發(fā)并測試了所提出的自動定位技術(shù)。每位病人使用4個標(biāo)記，共計(jì)168幅標(biāo)記圖像，每套CT圖像和MR圖像平均處理時(shí)間分別為87s(時(shí)間范圍82～91s)和18秒(9～48s)。

從24幅CT圖像中識別出23個候選基準(zhǔn)點(diǎn)，即算法對CT圖像的標(biāo)記誤判率為4.66%。在144幅MR圖像中選擇具有最高平均強(qiáng)度的四個標(biāo)記后，從144個真實(shí)標(biāo)記中識別出131個，算法對MR圖像的標(biāo)記誤判率為9%。

在含有45種MR模態(tài)和11種CT模態(tài)的18種案例下測試了本文的方法。測試中僅使用常規(guī)MR圖像(切面尺寸：256*256或 512*512像素)和CT圖像(切面尺寸：512*512像素)，所有的圖像數(shù)據(jù)表示為DICOM圖像序列，每個文件均為軸向掃描切面的圖像文件。經(jīng)測試，得出了理想結(jié)果。因文章篇幅限制，表1給出了兩種案例的具體結(jié)果。

表1 案例圖像信息和基準(zhǔn)點(diǎn)測試結(jié)果

使用CT或MR掃描儀獲取的圖像質(zhì)量的差異，并非所有模態(tài)下的基準(zhǔn)標(biāo)記均是可見的。例如，盡管MR的T1和T2模態(tài)下全部的12個基準(zhǔn)標(biāo)記均可見，但在案例1中的CT模態(tài)下，僅能通過視覺識別出9個基準(zhǔn)標(biāo)記。

并沒有對CT模態(tài)下漏掉的3個基準(zhǔn)標(biāo)記進(jìn)行重建，原因如下：1)認(rèn)為所漏掉的標(biāo)記在CT模態(tài)下并不存在，因?yàn)殛P(guān)于它們并沒有可用的信息，如像素強(qiáng)度等；2)因?yàn)榛鶞?zhǔn)標(biāo)記檢測的目的在于實(shí)現(xiàn)自動配置而非相反目的，不能使用其他模態(tài)下的基準(zhǔn)標(biāo)記位置信息對它們進(jìn)行CT模態(tài)重建。

表1中最后一行為實(shí)驗(yàn)結(jié)果。在兩個案例中，本文方法均檢測出了所有可見基準(zhǔn)標(biāo)記，沒有基準(zhǔn)標(biāo)記漏檢，檢測到的標(biāo)記位置均無誤。

在全部55個模態(tài)下(MR和CT)，共計(jì)430個基準(zhǔn)點(diǎn)(其中107個在CT圖像中)，僅在一個MR模態(tài)下的一個基準(zhǔn)點(diǎn)未能檢測出來，未檢測到的基準(zhǔn)點(diǎn)的軸位視圖如圖13所示，該基準(zhǔn)標(biāo)記并沒有位于圖像的頭骨部分，因此對此圖像序列的邊緣檢測沒有生成任何可識別模板。

圖13 本文方法未能檢測到的基準(zhǔn)標(biāo)記

4 結(jié)論

本文提出了一種全新的全自動基準(zhǔn)點(diǎn)檢測算法，可用于不同的圖像模態(tài)中，經(jīng)18個案例測試，實(shí)驗(yàn)結(jié)果準(zhǔn)確，誤判率很低，這個方法有助于實(shí)現(xiàn)基于基準(zhǔn)點(diǎn)的全自動圖像配準(zhǔn)。該方法和人眼識別具有相同的性能，不受像素強(qiáng)度值，基準(zhǔn)標(biāo)記位置和方位的影響。后續(xù)關(guān)于基準(zhǔn)點(diǎn)檢測的研究將關(guān)注以下方面：將基準(zhǔn)點(diǎn)檢測與圖像配準(zhǔn)向結(jié)合； 測試基準(zhǔn)點(diǎn)檢測方法對iMRI圖像的應(yīng)用。