顱內電極定位中的Zernike矩快速算法研究

劉思偉，董碩，嚴漢民

首都醫科大學宣武醫院 醫學工程處，北京 100053

顱內電極定位中的Zernike矩快速算法研究

劉思偉，董碩，嚴漢民

首都醫科大學宣武醫院 醫學工程處，北京 100053

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

幾何中心法；Hough變換法；多模圖像配準；顱內電極定位；Zernike矩

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

近年來，Zernike矩的方法在醫學圖像配準領域得到廣泛應用，其具有特征相關性小、冗余性小、抗噪能力強等特點[3-4]。通過改進有關矩的數學特性，Zernike矩的方法已經能夠達到亞像素級的配準效果[5]。

在實際應用Zernike矩的方法對圖像進行配準時，常遇到數據量巨大、計算時間過長等問題，且該矩要求待處理數據全部需在單位圓內。基于此，本研究通過應用幾何中心法和Hough變換法來提取圖像的特征，以獲取待配準的兩幅圖像各自的幾何中心坐標和配準區域，以縮小Zernike矩在計算最佳平移參數和伸縮比例時的計算范圍，尋找最佳單位圓范圍，從而減小計算耗時，提升配準效果。

1 研究方法

1.1 算法原理

幾何中心法是數學領域里描述形狀特征的方法。本研究將圖像特征區域的質心定義為幾何中心，以幾何中心為圓心所構成的圓形區域為配準區域。在計算時采用加權法和非加權法兩種方法來獲取質心，之后以該質心為圓心，計算特征半徑，以得到配準區域[6-7]。

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

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

圖1 模型圖像

1.2 處理方法

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

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

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

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

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

2 結果

2.1 實驗數據

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

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

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

2.2 測試結果

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

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

圖3 模型圖像幾何中心法測試結果

圖4 模型圖像Hough變換法測試結果

圖5 MR I圖像幾何中心法計算結果

圖6 CT圖像幾何中心法計算結果

圖7 CT和MR I圖像Hough變換法計算結果

3 討論

本研究以在包含顱腦全部信息的基礎上所得到的配準區域最小的幾何參數和配準區域為最優結果。對于質心相同或半徑相同的情況，以配準區域小或覆蓋噪聲區域小為更優結果。

由圖3和表1可知，幾何中心法對于同組圖像中雙圓環或三圓環的識別準確率，有灰度加權的結果優于沒有灰度加權的結果；對于同組圖像，圓環數量越少，計算結果越準確，且近似同心圓組的準確率優于非同心圓組；單圓環圖像的測試結果均較好，與是否有灰度加權無關。

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

表1 模型圖像幾何中心法測量結果

表2 模型圖像Hough變換法測量結果

表3 MR I圖像幾何中心法計算結果

表4 CT圖像幾何中心法計算結果

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

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

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

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

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

4 結論

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

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

[1] 柳淵.顱內電極定位和三維顯示方法的研究[D].北京:首都醫科大學,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矩的不變性與計算實現[J].內蒙古農業大學學報,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] 王炎玲,杜建軍,郭新宇.基于力矩主軸和互信息的黃瓜莖切片圖像配準[J].中國體視學與圖像分析,2013,(2):115-122.

[8] 羅述謙,周果宏.醫學圖像處理與分析[M].2版.北京:科學出版社,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] 盧蓉,范勇,陳念年.一種提取目標圖像最小外接矩形的快速算法[J].計算機工程,2010,36(21):178-180.

[12] 王琳婧,張書旭,林生趣,等.基于Viscous Fluid模型的快速CTCBCT圖像變形配準算法研究[J].中國醫療設備,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）；國家

自然科學基金項目資助（81372923）。

嚴漢民，教授，碩士生導師。

作者郵箱：liusw911@hotmail.com