林迪逵，吳劍威，李威，付麗媛，熊暉，陳自謙

解放軍福州總醫(yī)院 a. 醫(yī)學(xué)工程科/福建省醫(yī)學(xué)裝備管理質(zhì)量控制中心；b. 醫(yī)學(xué)影像中心，福建 福州 350025

引言

MRI掃描結(jié)果的精確與否是評(píng)判MRI技術(shù)臨床應(yīng)用成果的基礎(chǔ)標(biāo)準(zhǔn)，我們希望MRI圖像能夠?qū)颊叩慕馄?、病理和生理功能作出如?shí)呈現(xiàn)，但是問(wèn)題在于無(wú)論對(duì)于臨床研究或是診斷，實(shí)現(xiàn)高質(zhì)量圖像的獲取都是一項(xiàng)艱巨的任務(wù)?；谏鲜雠R床應(yīng)用背景，MRI研究與臨床工作者提出了質(zhì)量控制（Quality Control，QC）概念，旨在對(duì)磁共振成像質(zhì)量實(shí)現(xiàn)監(jiān)測(cè)，并通過(guò)一系列的技術(shù)程序保證MRI圖像質(zhì)量的持續(xù)提高。MRI質(zhì)量控制程序圍繞圖像質(zhì)量參數(shù)的檢測(cè)進(jìn)行展開(kāi)，通過(guò)成像質(zhì)量來(lái)客觀(guān)反映成像設(shè)備的運(yùn)行狀態(tài)，實(shí)現(xiàn)對(duì)成像系統(tǒng)性能狀態(tài)的評(píng)價(jià)。質(zhì)量控制程結(jié)果幫助臨床工程師對(duì)設(shè)備性能的改變作出反應(yīng)，保證MRI成像系統(tǒng)在安全、穩(wěn)定的狀態(tài)下運(yùn)行[1]。

MRI QC程序?qū)τ诳臻g線(xiàn)性度的檢測(cè)目的是為驗(yàn)證MRI系統(tǒng)是否以反映研究對(duì)象真實(shí)尺寸的縮放方式進(jìn)行成像。良好的空間線(xiàn)性度保證MRI圖像對(duì)成像目標(biāo)距離進(jìn)行真實(shí)再現(xiàn)，實(shí)現(xiàn)對(duì)組織幾何關(guān)系的精確反映，有助于臨床醫(yī)師進(jìn)行準(zhǔn)確的影像診斷。影響磁共振成像空間線(xiàn)性度的因素有很多，不合適的梯度縮放校正因子、不均勻磁場(chǎng)等成像環(huán)境都能導(dǎo)致MRI圖像空間線(xiàn)性度的下降[2-3]。

確定空間線(xiàn)性度的傳統(tǒng)方法依靠人工檢測(cè)，其存在以下缺陷：① 涉及多組重復(fù)測(cè)距，過(guò)程繁瑣耗時(shí)；② 檢測(cè)精度受限于QC檢測(cè)人員感覺(jué)器官的生理變化；③ 若檢測(cè)人員發(fā)生更替，QC的連續(xù)性意義也將大打折扣。為解決上述問(wèn)題，本文擬設(shè)計(jì)空間線(xiàn)性度檢測(cè)算法，取代傳統(tǒng)人工檢測(cè)，實(shí)現(xiàn)空間線(xiàn)性度的自動(dòng)測(cè)量，保證參數(shù)檢測(cè)結(jié)果如實(shí)反映成像設(shè)備相關(guān)性能狀態(tài)。

1 材料與方法

1.1 Magphan SMR 170體模的空間線(xiàn)性度檢測(cè)方法

Magphan SMR 170體模由美國(guó)體模實(shí)驗(yàn)室研制，該體模在橫斷位上設(shè)計(jì)5個(gè)測(cè)試層面，層面上設(shè)置的多組結(jié)構(gòu)模塊允許MRI QC程序?qū)举|(zhì)量控制參數(shù)進(jìn)行檢測(cè)，達(dá)到校驗(yàn)成像系統(tǒng)的目的。

QC程序?qū)臻g線(xiàn)性度的檢測(cè)在體模第四層橫斷位圖像上進(jìn)行（圖1）。該過(guò)程主要針對(duì)測(cè)試層面上若干圓孔間的間距進(jìn)行測(cè)量，定義空間線(xiàn)性度l為：

其中LR為體模中實(shí)際圓孔間距，三組孔間實(shí)際距離分別為80、40、20 mm；LM為人工檢測(cè)距離[4-5]。

圖1 基于Magphan SMR 170體模的空間線(xiàn)性度檢測(cè)示例

從上述檢測(cè)過(guò)程可以看出，若要實(shí)現(xiàn)空間線(xiàn)性度的自動(dòng)檢測(cè)，關(guān)鍵在于對(duì)圖像中若干目標(biāo)模塊（即圓孔）的識(shí)別，同時(shí)圓孔定位的精確性決定了線(xiàn)性度測(cè)量結(jié)果的準(zhǔn)確性。

1.2 基于隨機(jī)Hough變換的圓檢測(cè)

在數(shù)字圖像的自動(dòng)檢測(cè)過(guò)程中，對(duì)目標(biāo)形狀的檢測(cè)要求常常作為子問(wèn)題出現(xiàn)，例如實(shí)現(xiàn)對(duì)直線(xiàn)、圓或橢圓的識(shí)別。在大多數(shù)情況下，在預(yù)處理階段通過(guò)一些邊緣檢測(cè)程序的實(shí)施，可以獲得圖像空間中期望曲線(xiàn)上的圖像點(diǎn)或圖像像素。但是由于圖像數(shù)據(jù)質(zhì)量較差或者邊緣檢測(cè)程序存在缺陷等原因，可能存在期望曲線(xiàn)上的點(diǎn)或像素缺失等情況，導(dǎo)致上述邊緣噪聲點(diǎn)與理想的目標(biāo)形狀將存在空間偏差。所以，為達(dá)到目標(biāo)形狀的檢測(cè)要求，將提取到的邊緣特征聚類(lèi)成一系列合適的目標(biāo)形狀就顯得尤為重要，霍夫變換（Hough Transform）即為解決上述問(wèn)題而提出的。Hough變換針對(duì)參數(shù)化的圖像對(duì)象，通過(guò)設(shè)置累加器執(zhí)行“投票”過(guò)程，對(duì)峰值點(diǎn)作出判斷，從而實(shí)現(xiàn)將一系列邊緣點(diǎn)列為候選對(duì)象的目標(biāo)。

基于標(biāo)準(zhǔn)Hough變換（Standard Hough Transform，SHT）的圓檢測(cè)算法由于引進(jìn)了高維參數(shù)空間，在某些圖像實(shí)時(shí)處理的應(yīng)用環(huán)境中，變換效果并不理想。為解決SHT算法缺陷，隨機(jī)Hough變換（Randomized Hough Transform，RHT）應(yīng)運(yùn)而生。RHT的特點(diǎn)在于采用了隨機(jī)抽樣機(jī)制，實(shí)現(xiàn)圖像空間到參數(shù)空間的收斂映射，即隨機(jī)選取圖像空間上不共線(xiàn)的3個(gè)點(diǎn)，將上述三點(diǎn)映射至參數(shù)空間的一個(gè)點(diǎn)。相對(duì)于SHT的“一到多”映射，RHT提出的“多到一”映射無(wú)論對(duì)于內(nèi)存的占用，還是對(duì)于時(shí)間的消耗，都在很大程度上節(jié)省了開(kāi)銷(xiāo)[6-7]。RHT圓檢測(cè)原理如下。

隨機(jī)選取待識(shí)別圓上不共線(xiàn)的三點(diǎn)（x1，y1）（x2，y2）（x3,y3），分別代入圓方程可得方程組：

以上方程組能確定一個(gè)圓心為（a0，b0），半徑為r0的圓o0。下面的任務(wù)在于確定該待識(shí)別圓是否為真實(shí)的圓。取第四個(gè)圓上點(diǎn)（x4， y4），帶入圓心為（a0，b0）的圓方程，得到圓半徑1r′。將r0、1r′代入下式：

設(shè) ε0為預(yù)設(shè)誤差值。若 ε1＜ε0，認(rèn)為點(diǎn)（x4， y4）在圓 o0上，o0確定為候選圓。將待識(shí)別圓上的其余點(diǎn)（xi， yi）下的εi值，若εi＜ε0，累加值k加1，當(dāng)k達(dá)到預(yù)設(shè)閾值kthreshold時(shí)，則認(rèn)為o0為真實(shí)的圓。

上述RHT隨機(jī)采樣機(jī)制的一個(gè)特點(diǎn)在于采樣的無(wú)目標(biāo)性，雖然相較于SHT它在算法結(jié)構(gòu)上作出了有效的改進(jìn)，但是無(wú)法避免無(wú)效單元的引進(jìn)（圖2）。圖2中，o1、o2為待識(shí)別的目標(biāo)圓；d1、d2、d3為三個(gè)隨機(jī)抽樣點(diǎn)，其位于不同目標(biāo)圓上且不共線(xiàn)。根據(jù)RHT原理，認(rèn)為上述3個(gè)點(diǎn)可以確定候選圓o3，其圓參數(shù)為p，之后將進(jìn)行搜索參數(shù)單元集p等一系列程序來(lái)確定圓o3是否為真實(shí)圓。由于圓o3為虛假圓，所以后續(xù)對(duì)于圓o3的檢測(cè)程序均為無(wú)效累積過(guò)程。

圖2 RHT無(wú)效累積示例

假設(shè)有N個(gè)待識(shí)別圓，q為每個(gè)待識(shí)別圓上所選取點(diǎn)的數(shù)量，n為非圓上點(diǎn)。則隨機(jī)選擇的3點(diǎn)落在同一個(gè)待識(shí)別圓上的概率為：

如果剔除非圓點(diǎn)n，上式簡(jiǎn)化為：

本文空間線(xiàn)性度待檢測(cè)圖像上待識(shí)別的圓孔一共有13個(gè)，如果不考慮體模圖像中非目標(biāo)圓上的點(diǎn)，由式(7)可知，隨機(jī)選擇的3點(diǎn)落在同一個(gè)待識(shí)別圓孔上的概率為1/169，無(wú)效累積概率為168/169，如果將非圓邊緣點(diǎn)考慮在內(nèi)，無(wú)效累積現(xiàn)象將更為嚴(yán)重。針對(duì)本文圓檢測(cè)的應(yīng)用環(huán)境，認(rèn)為可以對(duì)上述RHT算法作出改進(jìn)，重點(diǎn)解決無(wú)效累積問(wèn)題，進(jìn)一步提高圓檢測(cè)效率。

1.3 改進(jìn)的RHT圓檢測(cè)算法

1.3.1 利用梯度方向信息減少無(wú)效累積

o1、o2為待識(shí)別的目標(biāo)圓；d1、d2、d3為三個(gè)隨機(jī)抽樣點(diǎn)，其位于不同目標(biāo)圓上且不共線(xiàn)（圖3）。由于d1與d2由于位于圓o1上，故d1與d2梯度方向所在的直線(xiàn)與必經(jīng)過(guò)圓心o1，且d1與d2連線(xiàn)的中垂線(xiàn)也必經(jīng)過(guò)圓心o1。反觀(guān)d1與d3，雖直線(xiàn)與將相交于點(diǎn)d4，但d1與d3連線(xiàn)的中垂線(xiàn)將不再經(jīng)過(guò)點(diǎn)d4，由此可以認(rèn)定d1與d3不在同一個(gè)圓?；谏鲜雠袛?，后續(xù)算法將不再進(jìn)行搜索參數(shù)單元集P等一系列參數(shù)空間累積過(guò)程來(lái)確定d1與d3是否同圓，省去大量的無(wú)效累積過(guò)程[8-9]。

圖3 利用梯度方向信息的改進(jìn)RHT示例

另一方面，由于改進(jìn)的RHT引入梯度方向信息進(jìn)行判斷，隨機(jī)抽樣點(diǎn)的個(gè)數(shù)從3個(gè)減少為2個(gè)，不考慮非圓點(diǎn)n，可知兩點(diǎn)落在同一個(gè)待識(shí)別圓上的概率為：

可以看出，相對(duì)于3點(diǎn)，保證2點(diǎn)共圓的概率大幅增加，對(duì)于存在大量待識(shí)別圓的檢測(cè)環(huán)境，這種概率的增加更為明顯，無(wú)效累積現(xiàn)象明顯減少。

1.3.2 利用先驗(yàn)知識(shí)減少無(wú)效累積

空間線(xiàn)性度檢測(cè)層面上除了待識(shí)別的圓孔，還包含1個(gè)體模輪廓圓以及4個(gè)支撐柱圓，為剔除無(wú)效信息對(duì)圓檢測(cè)過(guò)程的干擾，利用基于種子區(qū)域生長(zhǎng)法對(duì)體模中心感興趣區(qū)域進(jìn)行提取，實(shí)現(xiàn)感興趣區(qū)域的分割。分割過(guò)程如下：對(duì)圖像進(jìn)行二值化預(yù)處理后，在感興趣區(qū)域設(shè)置初始種子像素點(diǎn)（圖4）；通過(guò)檢查初始種子點(diǎn)相鄰像素的灰度特征，考慮是否將該相鄰像素添加至生長(zhǎng)區(qū)域，若種子與相鄰像素的灰度差大于設(shè)定閾值，區(qū)域停止生長(zhǎng)，得到體模中心感興趣區(qū)域（圖5）[10]。

由于感興趣區(qū)域中還包含條狀物模塊，圓檢測(cè)過(guò)程中也將會(huì)對(duì)其邊緣點(diǎn)集進(jìn)行提取。為避免累積無(wú)效圓參數(shù)，下面將利用半徑的判斷條件決定是否對(duì)圓參數(shù)p進(jìn)行保留。

d1(x1,y1)為圓邊緣點(diǎn)，d2(x2,y2)為矩形邊緣點(diǎn)，兩點(diǎn)梯度直線(xiàn)為l1、l2，兩點(diǎn)連線(xiàn)中垂線(xiàn)為l3，此時(shí)l1、l2、l3交于點(diǎn)c(xc, yc)，若只考慮梯度信息，則認(rèn)為兩點(diǎn)可能同圓，并將候選圓參數(shù)p加入p，造成誤判（圖5）。本文將引入半徑預(yù)設(shè)閾值o1，若邊緣點(diǎn)d1(x1,y1)、d2(x2,y2)與交點(diǎn)c(xc, yc)的距離的距離大于預(yù)設(shè)閾值o1，則認(rèn)為本次所選兩點(diǎn)不共圓，無(wú)需對(duì)該圓參數(shù)進(jìn)行累積，以提高檢測(cè)效率（圖6）。

圖4 設(shè)置初始種子像素點(diǎn)

圖5 種子生長(zhǎng)結(jié)束所分割的感興趣區(qū)域

圖6 梯度信息的誤判示例

綜合上述改進(jìn)措施，本文提出改進(jìn)的RHT（RHT+）圓檢測(cè)算法具體步驟，見(jiàn)圖7。

圖7 改進(jìn)的RHT（RHT+）圓檢測(cè)流程

1.4 空間線(xiàn)性度自動(dòng)檢測(cè)程序的實(shí)現(xiàn)

基于隨機(jī)Hough變換圓檢測(cè)算法的MRI QC空間線(xiàn)性度自動(dòng)檢測(cè)程序設(shè)計(jì)流程如下：

（1）對(duì)體模第四層橫斷位圖像進(jìn)行預(yù)處理。利用Otsu算法對(duì)圖像進(jìn)行二值化處理；針對(duì)圖像上出現(xiàn)的噪點(diǎn)，可采用數(shù)學(xué)形態(tài)學(xué)方法（即選取合適尺度的結(jié)構(gòu)元素，對(duì)圖像進(jìn)行膨脹、腐蝕操作）對(duì)其進(jìn)行消除[11]。

（2）提取圖像中模塊的邊緣信息。邊緣檢測(cè)算法利用邊緣增強(qiáng)算子對(duì)邊緣信息進(jìn)行強(qiáng)化，并對(duì)像素邊緣強(qiáng)度進(jìn)行計(jì)算，最后通過(guò)與閾值的比較，實(shí)現(xiàn)對(duì)邊緣點(diǎn)集的獲取。本文采用邊緣檢測(cè)Canny算子進(jìn)行邊緣檢測(cè)：應(yīng)用一階差分卷積模板實(shí)現(xiàn)目標(biāo)邊緣的增強(qiáng)；應(yīng)用非極大值抑制實(shí)現(xiàn)局部梯度方向最大值的保留；應(yīng)用雙線(xiàn)性閾值實(shí)現(xiàn)邊緣的檢測(cè)和連接[12]。

（3）分別利用RHT與RHT+圓檢測(cè)算法對(duì)目標(biāo)圓孔進(jìn)行檢測(cè)，并對(duì)識(shí)別到的圓孔進(jìn)行標(biāo)記。

（4）已知目標(biāo)圓孔位置，根據(jù)圖片尺寸度量由式(1)分別計(jì)算頻率與相位編碼方向上的空間線(xiàn)性度。

2 結(jié)果

選擇某次質(zhì)量控制得到的Magphan SMR 170體模第二層圖像作為檢測(cè)對(duì)象，圖像大小為256×256。算法編程環(huán)境為：MATLAB 2012b；運(yùn)行硬件環(huán)境為：CPU 2.3 GHz、內(nèi)存2 GB。

RHT與RHT+圓檢測(cè)算法對(duì)目標(biāo)圓孔的識(shí)別結(jié)果（識(shí)別到的圓孔用紅圈標(biāo)出），見(jiàn)圖8～9；RHT+與RHT圓檢測(cè)效果進(jìn)行比較，見(jiàn)表1；人工（由熟練的質(zhì)量控制工程師執(zhí)行此次空間線(xiàn)性度檢測(cè)程序）以及基于RHT+圓檢測(cè)算法的線(xiàn)性度測(cè)量結(jié)果，見(jiàn)表2。

圖8 RHT算法圓檢測(cè)效果

圖9 RHT+算法圓檢測(cè)效果

表1 RHT+與RHT圓檢測(cè)效果比較

綜合上述實(shí)驗(yàn)結(jié)果：RHT與RHT+圓檢測(cè)算法都實(shí)現(xiàn)了對(duì)13個(gè)目標(biāo)圓孔的識(shí)別，相較于RHT算法，RHT+由于加入剔除無(wú)效累積的改進(jìn)手段，檢測(cè)速度獲得近5倍的提升。表2中，基于RHT+算法的自動(dòng)檢測(cè)程序給出了空間線(xiàn)性度的測(cè)量結(jié)果，但該結(jié)果與人工檢測(cè)結(jié)果存在一定出入，為比較二者差異，本文采用統(tǒng)計(jì)學(xué)方法行差異顯著性檢驗(yàn)。這里提出H0假設(shè)：自動(dòng)檢測(cè)結(jié)果與人工檢測(cè)沒(méi)有差異，顯著性水平α設(shè)置為0.05級(jí)的情況下，P＞0.05時(shí)接受原假設(shè)，P＜0.05時(shí)拒絕原假設(shè)；檢驗(yàn)結(jié)果得到P=0.10009＞0.05，接受H0假設(shè)，即認(rèn)為兩種方法檢測(cè)結(jié)果差異無(wú)統(tǒng)計(jì)學(xué)意義[13]。

3 討論

近幾年，我國(guó)的MRI質(zhì)量控制相關(guān)工作與研究呈現(xiàn)蓬勃發(fā)展之勢(shì)，當(dāng)中對(duì)于磁共振質(zhì)量控制參數(shù)自動(dòng)檢測(cè)的討論也愈發(fā)深入。作者官能成等[14]2015年在雜志發(fā)表質(zhì)量控制自動(dòng)檢測(cè)相關(guān)研究成果，可以說(shuō)是開(kāi)創(chuàng)該研究領(lǐng)域之先河；直至2017年，張春青等[15]、李成等[16]重審研究背景，各自闡述了對(duì)質(zhì)量控制自動(dòng)檢測(cè)系統(tǒng)構(gòu)建方法的理解。以上等人的研究方向主要著眼于系統(tǒng)框架設(shè)計(jì)，未對(duì)參數(shù)檢測(cè)所涉及的圖像處理算法做具體討論。而針對(duì)空間線(xiàn)性度，李馥然[17]利用邊緣檢測(cè)以及連通域標(biāo)記等圖像處理方法來(lái)實(shí)現(xiàn)該參數(shù)的自動(dòng)檢測(cè)，孟奧等[18]則利用人工手動(dòng)選取目標(biāo)區(qū)域法對(duì)線(xiàn)性度檢測(cè)模塊坐標(biāo)進(jìn)行定位。相較于以上研究方法，本文針對(duì)空間線(xiàn)性度提出了具體的檢測(cè)算法，著力于檢測(cè)過(guò)程人工的完全解放，在保證參數(shù)檢測(cè)準(zhǔn)確度的同時(shí)，強(qiáng)調(diào)檢測(cè)效率，這對(duì)于后期構(gòu)建全參數(shù)、實(shí)時(shí)的磁共振質(zhì)量控制自動(dòng)檢測(cè)系統(tǒng)具有重要意義。

前文結(jié)果中所行差異顯著性檢驗(yàn)表明，兩種方法檢測(cè)結(jié)果差異無(wú)統(tǒng)計(jì)學(xué)意義，可以認(rèn)為本文所設(shè)計(jì)的自動(dòng)檢測(cè)程序能夠取代人工檢測(cè)方式進(jìn)行空間線(xiàn)性度的檢測(cè)，這對(duì)于MRI質(zhì)量控制程序的實(shí)施具備以下意義：① 摒棄了人工過(guò)程干預(yù)，無(wú)須受限于質(zhì)量控制檢測(cè)人員的檢測(cè)行為習(xí)慣、操作熟練程度等主觀(guān)因素影響，在追求隨機(jī)誤差最小化的過(guò)程中具備天然優(yōu)勢(shì)；② 算法運(yùn)行的穩(wěn)定性與連續(xù)性滿(mǎn)足了質(zhì)量控制程序?qū)τ跈z測(cè)連續(xù)性的要求，避免因?yàn)闄z測(cè)人員發(fā)生更替等原因而引起連續(xù)性的破壞，保證參數(shù)檢測(cè)結(jié)果客觀(guān)反映MRI系統(tǒng)相關(guān)性能狀態(tài)的真實(shí)波動(dòng)，這對(duì)于發(fā)現(xiàn)、隔離、解決成像系統(tǒng)細(xì)微問(wèn)題具有重大意義；③ 規(guī)避質(zhì)控程序?qū)筇幚砉ぷ髡鹃L(zhǎng)時(shí)間占用的現(xiàn)象，解決系統(tǒng)校驗(yàn)行為與臨床使用需求的沖突。

表2 人工與基于RHT+圓檢測(cè)算法的線(xiàn)性度測(cè)量結(jié)果（%）

對(duì)于算法結(jié)構(gòu)，本文所提出的空間線(xiàn)性度自動(dòng)檢測(cè)方法還應(yīng)用了包括圖像預(yù)處理、圖像增強(qiáng)等算法，由于篇幅原因，本文對(duì)上述圖像處理技術(shù)未做詳細(xì)討論，但不能否認(rèn)其對(duì)自動(dòng)檢測(cè)目標(biāo)實(shí)現(xiàn)的重要作用。此外，本文研究雖取得一定的成果，但若在更為復(fù)雜、嚴(yán)苛際成像環(huán)境中，如何在保證參數(shù)檢測(cè)準(zhǔn)確性的同時(shí)，增強(qiáng)算法的穩(wěn)定性、魯棒性、快速性，這是本文后續(xù)的研究目標(biāo)，也是其他質(zhì)量控制相關(guān)參數(shù)要實(shí)現(xiàn)自動(dòng)檢測(cè)的設(shè)計(jì)理念。