劉思偉，董碩，嚴(yán)漢民

首都醫(yī)科大學(xué)宣武醫(yī)院 醫(yī)學(xué)工程處，北京 100053

顱內(nèi)電極定位中的Zernike矩快速算法研究

劉思偉，董碩，嚴(yán)漢民

首都醫(yī)科大學(xué)宣武醫(yī)院 醫(yī)學(xué)工程處，北京 100053

目的 減少顱內(nèi)電極定位中多模圖像Zernike矩的配準(zhǔn)耗時，提升其配準(zhǔn)精度。方法 用幾何中心法和Hough變換法對人腦CT和MRI圖像進行幾何參數(shù)和配準(zhǔn)區(qū)域的計算，并根據(jù)人腦結(jié)構(gòu)特點建立模型對算法有效性進行測試。本文分別應(yīng)用骨架算法、顱骨提取法、腐蝕聯(lián)合大津閾值法、邊緣增強法、邊緣增強聯(lián)合腐蝕法、Bottom hat法、Sobel邊緣提取法對人腦CT和MRI圖像進行預(yù)處理。結(jié)果 對于模型圖像，Hough變換法和幾何中心法測得的質(zhì)心與實際質(zhì)心坐標(biāo)分別相差±1個像素和±5個像素，主軸誤差范圍均＜3%。對于CT和MRI圖像，幾何中心法的計算結(jié)果更優(yōu)，其質(zhì)心點和主軸的計算誤差范圍分別可控制在周圍5鄰域內(nèi)和±6%內(nèi)。結(jié)論 幾何中心法和Hough變換法均能夠很好地計算出模型圖像幾何參數(shù)和配準(zhǔn)區(qū)域。對于CT和MRI圖像，幾何中心法的計算效果較好，Hough變換法的計算效果有待提升。

幾何中心法；Hough變換法；多模圖像配準(zhǔn)；顱內(nèi)電極定位；Zernike矩

癲癇是由多種病因引起的慢性腦部疾患，每年新發(fā)病人約30萬。其中約20%的病人不能采用藥物控制，而在這些病人之中，約50%的人可以采用手術(shù)治療，而致癇灶的準(zhǔn)確定位是手術(shù)成功的關(guān)鍵[1]。顱內(nèi)電極腦電記錄法具有靈敏度高、定位精確的優(yōu)點，可以對致癇灶進行準(zhǔn)確定位。但含有顱內(nèi)電極的CT圖像由于存在金屬偽影和金屬電極信息，因此與MRI圖像的配準(zhǔn)效果欠佳，最終影響致癇灶的定位[2]。

近年來，Zernike矩的方法在醫(yī)學(xué)圖像配準(zhǔn)領(lǐng)域得到廣泛應(yīng)用，其具有特征相關(guān)性小、冗余性小、抗噪能力強等特點[3-4]。通過改進有關(guān)矩的數(shù)學(xué)特性，Zernike矩的方法已經(jīng)能夠達(dá)到亞像素級的配準(zhǔn)效果[5]。

在實際應(yīng)用Zernike矩的方法對圖像進行配準(zhǔn)時，常遇到數(shù)據(jù)量巨大、計算時間過長等問題，且該矩要求待處理數(shù)據(jù)全部需在單位圓內(nèi)?；诖耍狙芯客ㄟ^應(yīng)用幾何中心法和Hough變換法來提取圖像的特征，以獲取待配準(zhǔn)的兩幅圖像各自的幾何中心坐標(biāo)和配準(zhǔn)區(qū)域，以縮小Zernike矩在計算最佳平移參數(shù)和伸縮比例時的計算范圍，尋找最佳單位圓范圍，從而減小計算耗時，提升配準(zhǔn)效果。

1 研究方法

1.1 算法原理

幾何中心法是數(shù)學(xué)領(lǐng)域里描述形狀特征的方法。本研究將圖像特征區(qū)域的質(zhì)心定義為幾何中心，以幾何中心為圓心所構(gòu)成的圓形區(qū)域為配準(zhǔn)區(qū)域。在計算時采用加權(quán)法和非加權(quán)法兩種方法來獲取質(zhì)心，之后以該質(zhì)心為圓心，計算特征半徑，以得到配準(zhǔn)區(qū)域[6-7]。

Hough變換法利用圖像的全局特征直接檢測目標(biāo)輪廓，在預(yù)先了解區(qū)域大致形狀的基礎(chǔ)上，在圖像中尋找預(yù)期的形狀，因此可以作為簡單的模板匹配技術(shù)對圖像對象進行識別[8]。Hough變換的主要優(yōu)點是受噪聲影響較小，測得結(jié)果為依概率分布的函數(shù)。本文選取的兩幅待配準(zhǔn)圖像均為人腦圖像，而人腦顱骨可以粗略認(rèn)為其近似圓形，符合Hough變換的適用條件[9]。

為驗證上述算法的有效性，本文建立了如下模型進行測試：分別采用單個、雙個和3個正圓圓環(huán)組合成的近似同心圓和非近似同心圓的兩組圖像模擬人腦（圖1）。人腦圖像可以近似成是由顱骨和腦組織這兩個質(zhì)心近似相同的圓環(huán)組成的，但對于腦部放有電極的患者而言，由于電極分布無序，且其與顱骨在CT圖像上的灰度值基本相近，因此在用灰度信息計算質(zhì)心時可能會造成顱骨與腦組織的質(zhì)心不再重合，即顱骨和腦組織變?yōu)閳A心不同的圓環(huán)，故本研究設(shè)計了上述兩組圖像對算法進行測試。與此同時，為了更好地模擬人腦，本研究對兩組圖像還做了灰度值設(shè)定。

圖1 模型圖像

1.2 處理方法

由于Hough變換計算耗時較大，因此本研究使用該方法測試模型圖像時，先利用邊緣檢測算子對圖像進行邊緣提取，再根據(jù)幾何中心法得到的結(jié)果設(shè)定Hough變換計算范圍，以減少計算耗時[10]。

考慮到圖像噪聲對算法的影響，本研究在應(yīng)用幾何中心法和Hough變換法對CT和MRI圖像進行計算前對相關(guān)圖像進行了預(yù)處理。使用幾何中心法對MRI圖像進行計算的步驟如下：① 采用大津閾值對原始圖像進行分割；② 使用腐蝕算法對分割結(jié)果進行腐蝕；③ 分別利用有無灰度加權(quán)的方法計算圖像質(zhì)心；④ 根據(jù)相鄰的4鄰域?qū)D像進行內(nèi)部點消除；⑤ 以幾何中心點為圓心計算圖像配準(zhǔn)區(qū)域特征半徑。使用幾何中心法對CT圖像進行計算的步驟如下：① 采用大津閾值對原始圖像進行分割；② 使用腐蝕算法對分割結(jié)果進行腐蝕；③ 分別采用有無灰度加權(quán)的方法計算圖像質(zhì)心；④ 對步驟②得到的圖像繼續(xù)進行腐蝕、內(nèi)部點消除和孤立點消除；⑤ 分別采用Bottom hat、邊緣增強、邊緣增強聯(lián)合腐蝕的方法處理圖像；⑥ 以幾何中心點為圓心計算圖像配準(zhǔn)區(qū)域特征半徑。

但在預(yù)實驗中，上述方法對圖像質(zhì)心的計算效果并不理想，因此對其進行了如下改進：① 提取原始圖像中的顱骨和金屬電極；② 使用剪影法去除原始圖像中的顱骨和金屬電極信息；③ 利用大津閾值對剪影圖像進行分割并計算幾何中心；④ 分別采用Bottom hat、邊緣增強、邊緣增強聯(lián)合腐蝕的方法計算半徑。

本研究在應(yīng)用Hough變換法計算CT與MRI圖像的幾何特征時，利用幾何中心法所得的質(zhì)心為圓心、所得半徑的不同倍數(shù)為半徑，構(gòu)建同心圓環(huán)，兩圓之間的區(qū)域為計算使用區(qū)域[11]，此設(shè)計可以在保證計算準(zhǔn)確性的基礎(chǔ)上減小計算量。對于Hough變換使用的CT和MRI數(shù)據(jù)，本研究在幾何中心法處理的基礎(chǔ)上分別進行了以下4種處理：① 提取圖像的顱骨部分作為輸入數(shù)據(jù)；② 進行腐蝕和大津閾值處理后作為輸入數(shù)據(jù)；③ 用Sobel邊緣提取法處理后作為輸入數(shù)據(jù)；④ 用骨架算法處理后作為輸入數(shù)據(jù)。

考慮到CT和MRI圖像中的偽影和部分噪聲的灰度值與顱骨的灰度值相同或相近[12]，且在空間上相互連接，因此本研究在提取顱骨信息時也會涵蓋上述信息。

2 結(jié)果

2.1 實驗數(shù)據(jù)

本研究的原始數(shù)據(jù)包括頭顱CT掃描圖像和頭顱MRI掃描圖像。其中頭顱CT圖像數(shù)據(jù)為：植入電極后的CT掃描斷層圖像，共227層，單幅圖像512像素×512像素，像素大小0.5 mm×0.5 mm，層厚1 mm。頭顱MRI掃描圖像數(shù)據(jù)為：植入電極前的軸位T1W1像，共268層，單幅圖像512像素×512像素，像素大小0.5 mm×0.5 mm，層厚1 mm。原始圖像均為DICOM格式數(shù)據(jù)集，由首都醫(yī)科大學(xué)宣武醫(yī)院提供。圖2為原始圖像示例。

圖2 CT、MR I原始圖像（左圖為植入電極后的CT圖像，右圖為植入電極前的MR I圖像）

本研究將模型圖像分為兩組，每組又分為各個圓環(huán)灰度值均相同的和均不同的兩組數(shù)據(jù)，灰度值各異組的灰度值從外到內(nèi)設(shè)定為64、32、16，兩組模型圖像的最外層圓環(huán)的圓心坐標(biāo)和半徑均分別為[194,234]和156 Pixel。

2.2 測試結(jié)果

本研究利用MATLAB2007a進行程序編寫和結(jié)果測試。Hough變換法以幾何中心法所得的圓心為圓心、所得半徑的1.02倍為半徑，構(gòu)建計算區(qū)域；并在原圖上繪制出經(jīng)計算得到的配準(zhǔn)區(qū)域和質(zhì)心，以便對兩種算法的計算結(jié)果進行比較。

圖3和表1為模型圖像幾何中心法的測試結(jié)果。圖4和表2為模型圖像Hough變換法的測試結(jié)果。圖5和表3為應(yīng)用幾何中心法對MRI圖像進行計算的結(jié)果。圖6和表4為應(yīng)用幾何中心法對CT圖像進行計算的結(jié)果。圖7和表5為使用Hough變換法對CT和MRI圖像進行計算的結(jié)果。

圖3 模型圖像幾何中心法測試結(jié)果

圖4 模型圖像Hough變換法測試結(jié)果

圖5 MR I圖像幾何中心法計算結(jié)果

圖6 CT圖像幾何中心法計算結(jié)果

圖7 CT和MR I圖像Hough變換法計算結(jié)果

3 討論

本研究以在包含顱腦全部信息的基礎(chǔ)上所得到的配準(zhǔn)區(qū)域最小的幾何參數(shù)和配準(zhǔn)區(qū)域為最優(yōu)結(jié)果。對于質(zhì)心相同或半徑相同的情況，以配準(zhǔn)區(qū)域小或覆蓋噪聲區(qū)域小為更優(yōu)結(jié)果。

由圖3和表1可知，幾何中心法對于同組圖像中雙圓環(huán)或三圓環(huán)的識別準(zhǔn)確率，有灰度加權(quán)的結(jié)果優(yōu)于沒有灰度加權(quán)的結(jié)果；對于同組圖像，圓環(huán)數(shù)量越少，計算結(jié)果越準(zhǔn)確，且近似同心圓組的準(zhǔn)確率優(yōu)于非同心圓組；單圓環(huán)圖像的測試結(jié)果均較好，與是否有灰度加權(quán)無關(guān)。

由圖4和表2可知，Hough變換法對于近似同心圓組和非同心圓組均有較好的識別準(zhǔn)確率，且不受圓環(huán)數(shù)量的影響；對于同一組圓環(huán)而言，圓環(huán)數(shù)量的差異幾乎不影響計算結(jié)果。因此，該方法對于模型圖像的測試結(jié)果要好于幾何中心法。與此同時，由于Hough變換法在計算中并不涉及灰度信息，因此不需要根據(jù)灰度信息進行分組測試。由圖4還可知，Hough變換法得到的半徑并非為最外層圓環(huán)的外半徑，而是其內(nèi)半徑。因此對于顱腦圖像而言，該方法可能無法得到顱腦最外層組織結(jié)構(gòu)所對應(yīng)的圓環(huán)半徑，這點在實際應(yīng)用時應(yīng)注意。由于Hough變換法識別圓形的計算結(jié)果實際上是圓心和半徑的依概率分布，對于以矩陣形式存儲的圓形，由于其邊緣無法絕對光滑，所以計算出的概率最大值可能并非連續(xù)域理論上的概率最大值。

表1 模型圖像幾何中心法測量結(jié)果

表2 模型圖像Hough變換法測量結(jié)果

表3 MR I圖像幾何中心法計算結(jié)果

表4 CT圖像幾何中心法計算結(jié)果

表5 CT和MR I圖像Hough變換法計算結(jié)果

由圖5和表3的幾何中心法計算結(jié)果可以看出，對于MRI圖像，加入腐蝕的結(jié)果好于未加入腐蝕的結(jié)果，未使用灰度加權(quán)的結(jié)果好于使用灰度加權(quán)的結(jié)果，帶有腐蝕且未加入灰度加權(quán)的方法可以較好地計算出MRI中顱腦的幾何中心和配準(zhǔn)區(qū)域。加入腐蝕的算法降低了顱骨外脂肪組織對計算的影響，因而其計算效果提升。而就灰度加權(quán)方法對結(jié)果的影響來看，模型圖像測試結(jié)果與真實圖像計算結(jié)果得出的結(jié)論相反，出現(xiàn)這種結(jié)果的原因是人腦顱骨并非為標(biāo)準(zhǔn)正圓且在CT圖像中的灰度值很大，因此較大的權(quán)重值和自身形狀反而成了影響測試結(jié)果的因素。

由圖6和表4的幾何中心法計算結(jié)果可以看出，對于CT圖像，在未去除金屬電極影響的情況下，測得的幾何中心應(yīng)偏靠金屬電極集中的位置。該方法計算的配準(zhǔn)區(qū)域和理想范圍相比偏靠電極集中且噪聲灰度值高的方向，因此會造成相反方向區(qū)域信息的丟失，其中采用邊緣增強聯(lián)合腐蝕方法處理的結(jié)果該現(xiàn)象最為明顯。3種方法計算出的配準(zhǔn)區(qū)域基本準(zhǔn)確，Bottom hat的效果最好，除右下角顱骨信息部分丟失外，保留了全部顱腦信息且引入的偽影和噪聲信息較少。

在去除金屬電極和顱骨信息之后，幾何中心法可以較好地計算出CT圖像的幾何中心，但是由于去除了顱骨信息，所以在半徑計算上出現(xiàn)了較大偏差。3種半徑計算方法中只有邊緣增強聯(lián)合腐蝕的方法能夠較好地計算出配準(zhǔn)區(qū)域，而另外兩種算法的計算誤差比較大。為此，本研究嘗試分別用去除顱骨和金屬電極信息的方法計算幾何中心點，用Bottom hat、邊緣增強、邊緣增強聯(lián)合腐蝕的方法計算未去除顱骨信息的圖像半徑，即聯(lián)合應(yīng)用兩種方法。從處理結(jié)果可以看出，在聯(lián)合應(yīng)用的3種方法里，邊緣增強方法計算的顱腦配準(zhǔn)區(qū)域除金屬電極集中方向上有信息少量缺失外，能夠包括其他全部顱腦信息。雖然聯(lián)合應(yīng)用方法比單純使用Bottom hat要更多地引入偽影和噪聲信息，但是其在保留顱腦信息方面的優(yōu)勢更為突出，能夠減少有效信息的損失，綜合處理效果更佳。

由圖7和表5可以看出，使用Hough變換法計算CT圖像與MRI圖像的幾何中心和半徑的效果較差，計算得到的幾何中心較理想值有所偏移，半徑較理想值偏小。對于CT圖像而言，Sobel邊緣提取法的計算結(jié)果最好，但是在顱骨邊緣區(qū)域仍有信息缺失，而其他3種方法得到的配準(zhǔn)區(qū)域顱腦信息缺失更多，相比于幾何中心法的最佳提取結(jié)果仍有不足；對于MRI圖像而言，Sobel邊緣提取法和顱骨提取法的效果相對較好，但是得到的配準(zhǔn)區(qū)域顱腦信息缺失嚴(yán)重，顱骨信息集中區(qū)域的顱骨信息大量缺失，并且有小部分腦組織信息缺失，而另外兩種方法得到的結(jié)果更差，與幾何中心法的最佳結(jié)果相比有很大差距。分析原因大致有以下幾點：① 所使用的Hough程序是按照標(biāo)準(zhǔn)正圓編寫的，而實際使用的圖像并非為標(biāo)準(zhǔn)正圓，因而會導(dǎo)致識別精度下降；② Hough變換法的計算結(jié)果為依概率的函數(shù)分布，因此最佳結(jié)果可能并非對應(yīng)概率最大時的函數(shù)值；③ 從模型測試圖像中可以看出Hough變換法在計算半徑時會得到內(nèi)側(cè)圓的半徑值，而不是外側(cè)圓的半徑值，因此在CT、MRI圖像計算中會出現(xiàn)半徑值較小的情況。

4 結(jié)論

相比于幾何中心法而言，Hough變換法的處理時間較長。對于模型圖像而言，這兩種方法均能取得較好的效果，并且Hough變換法不受灰度變化和干擾信息影響，在近似同心圓組和非近似同心圓組的圖像測試中均能準(zhǔn)確識別出最大圓的圓心，且計算結(jié)果與真實結(jié)果相同；但是在實際測試中，Hough變換法對于CT、MRI圖像的識別準(zhǔn)確率均不高，計算結(jié)果與真實值相差較大，但是幾何中心法則能很好地識別出顱腦區(qū)域。綜合而言，幾何中心法的計算結(jié)果優(yōu)于Hough變換法，其質(zhì)心點和主軸的計算范圍分別能控制在周圍5鄰域內(nèi)和±6%以內(nèi)。

與Hough變換法相比，幾何中心法對于CT、MRI圖像的識別準(zhǔn)確率較高，但是考慮到Hough變換法在識別時不受噪聲和偽影信息的影響，因而可以對Hough變換法進行改進以提升實際測量精度。除上述實驗外，對于CT圖像，能否利用金屬電極和電極偽影的不對稱性進行圖像幾何特征計算；對于CT、MRI圖像，不同的邊緣提取方法會如何影響Hough變換法的結(jié)果；對于矩和主軸的配準(zhǔn)方法，在獲得兩幅圖像幾何特征之后如何減小金屬電極和電極偽影對配準(zhǔn)效果的影響；除空間灰度外，引入時頻分析能否得到更好的配準(zhǔn)參數(shù)，對于這些問題我們也將展開深入研究，以更好地解決顱內(nèi)電極定位中多模圖像配準(zhǔn)問題[13]。

[1] 柳淵.顱內(nèi)電極定位和三維顯示方法的研究[D].北京:首都醫(yī)科大學(xué),2009.

[2] Miller KJ,Makeig S,Hebb AO,et al.Cortical electrode localization from X-rays and simple mapping for electrocorticographic research:the"Location on Cortex"(LOC) package for MATLAB[J]. Journal of Neuroscience Methods,2007,162(1-2):303-308.

[3] Yang Z,Fang T.On the accuracy of image normalization by Zernike moments[J].Image and Vision Computing,2010,28(3):403-413.

[4] 郝敏,麻碩士,侯振杰.Zernike矩的不變性與計算實現(xiàn)[J].內(nèi)蒙古農(nóng)業(yè)大學(xué)學(xué)報,2009,30(2):222-225.

[5] Ma Z,Kang B,Lv K,et al.Nonlinear radon transform using Zernike moment for shape analysis[J].Computational and Mathematical Methods in Medicine,2013,2013:1-9.

[6] Sotiras A,Davatzikos C,Paragios N.Deformable medical image registration:a survey[J].IEEE Transactions on Medical Imaging,2013,32(7):1153-1190.

[7] 王炎玲,杜建軍,郭新宇.基于力矩主軸和互信息的黃瓜莖切片圖像配準(zhǔn)[J].中國體視學(xué)與圖像分析,2013,(2):115-122.

[8] 羅述謙,周果宏.醫(yī)學(xué)圖像處理與分析[M].2版.北京:科學(xué)出版社,2010.

[9] Hermann E,Bleicken S,Subburaj Y,et al.Automated analysis of giant unilamellar vesicles using circular Hough transformation[J]. Bioinformatics,2014,30(12):1747-1754.

[10] Pauwels R,Jacobs R,Bosmans H,et al.Automated implant segmentation in cone-beam CT using edge detection and particle counting[J].International Journal of Computer Assisted Radiology and Surgery,2014,9(4):733-743.

[11] 盧蓉,范勇,陳念年.一種提取目標(biāo)圖像最小外接矩形的快速算法[J].計算機工程,2010,36(21):178-180.

[12] 王琳婧,張書旭,林生趣,等.基于Viscous Fluid模型的快速CTCBCT圖像變形配準(zhǔn)算法研究[J].中國醫(yī)療設(shè)備,2013,28(6): 21-23,33.

[13] Naranjo V,Lloréns R,Alca?iz M,et al.Metal artifact reduction in dental CT images using polar mathematical morphology[J].Computer Methods and Programs in Biomedicine,2011,102(1):64-74.

Research on Fast Algorithm for Zernike Moment in Intracranial Electrode Localization

LIU Si-wei, DONG Shuo, YAN Han-min
Department of Medical Engineering, Xuanwu Hospital, Capital Medical University, Beijing 100053, China

Objective To reduce the registration time and improve the registration accuracy of multimodality images with Zernike moment in intracranial electrode localization. Methods The geometric center method and Hough transform method were used in calculating the geometric parameters and registration areas of CT and MRI images of human brain. The validity of the two methods was tested with simulated models which were constructed according to the structural characteristics of human brain. Various methods including skeleton algorithm, skull extraction method, erode combined with OTSU threshold value method, edge enhancement method, edge enhancement combined with erode method, Bottom hat method and Sobel edge extraction method were used in the preprocessing of CT and MRI images. Results The center-of-mass coordinates of stimulated images measured by geometric center method and Hough transform method are different from real center-of-mass coordinates within ±1 pixel and ±5 pixels, respectively. The principle axis errors of the two methods are less than 3%. The results of relevant parameters of CT and MRI images calculated by geometric center method are better that those of CT and MRI images calculated by Hough transform method. The calculation errors of center-of-mass points and principle axis of CT and MRI images calculated by geometric center method can be controlled within 5 neighborhoods and ±6%, respectively. Conclusion Both geometric center method and Hough transform method can calculate the geometric parameters and registration areas of simulated models. However, the calculation effect of geometric center method is better than that of Hough transform method for the geometric parameters and registration areas of CT and MRI images.

geometric center method; Hough transform method; multimodality image registration; intracranial electrode localization; Zernike moment

R742.1；TH774

A

10.3969/j.issn.1674-1633.2015.01.006

1674-1633(2015)01-0023-05

2014-10-09

國家科技支撐計劃項目資助（2011BAI02B02）；國家

自然科學(xué)基金項目資助（81372923）。

嚴(yán)漢民，教授，碩士生導(dǎo)師。

作者郵箱：liusw911@hotmail.com