運(yùn)用投影正弦圖和點(diǎn)擴(kuò)展函數(shù)法確定CT系統(tǒng)旋轉(zhuǎn)中心＊

羅福婷，鄭繼明，張 亞，謝婷婷

（重慶郵電大學(xué) 理學(xué)院，重慶 400065）

計算機(jī)斷層成像（ComputedTomography，CT）在醫(yī)療檢測、無損檢測和逆向工程中發(fā)揮著重大的作用，被視為20世紀(jì)影響人類發(fā)展的重大技術(shù)之一。安裝CT系統(tǒng)時往往存在誤差，會影響成像質(zhì)量，因此，需要對安裝好的CT系統(tǒng)進(jìn)行參數(shù)標(biāo)定，根據(jù)標(biāo)定模板來確定CT系統(tǒng)旋轉(zhuǎn)中心是首要問題。本文基于2017年“高教社杯”全國大學(xué)生數(shù)學(xué)建模競賽A題[1]，討論如何精確計算一種典型的二維CT系統(tǒng)旋轉(zhuǎn)中心的問題。

CT系統(tǒng)旋轉(zhuǎn)中心定位的方法有幾何校正法、迭代校正法和投影正弦圖重心法等[2]。幾何校正法需已知精準(zhǔn)的中心射線束的位置、焦距，這在工程應(yīng)用中是不容易做到的，且會受計算效率的制約，實(shí)用性差；迭代校正法將旋轉(zhuǎn)中心偏移轉(zhuǎn)換為焦點(diǎn)問題，與手動調(diào)整原理相同，煩瑣、精度低。張蔚等[3]根據(jù)CT系統(tǒng)投影正弦圖的對稱性、正弦性、參數(shù)信息等，提出投影圖邊緣法，先對投影正弦圖進(jìn)行分割，并二值化投影正弦圖，將投影圖分割為空氣投影和工件投影兩部分，求物體投影在各個方向的質(zhì)心，最后將所有質(zhì)心位置加權(quán)求均值。該方法能達(dá)到10-1像素級，是理想的旋轉(zhuǎn)中心初始定位方法。為了進(jìn)一步提高CT系統(tǒng)旋轉(zhuǎn)中心定位的精度，本文提出了點(diǎn)擴(kuò)展函數(shù)法。

考慮到點(diǎn)擴(kuò)展函數(shù)僅能精確在1像素內(nèi)的偏移量，首先分析了模板投影正弦圖的特征，并對旋轉(zhuǎn)中心進(jìn)行初步定位，然后利用圖像重建確定點(diǎn)擴(kuò)展函數(shù)，從而較精確地確定出系統(tǒng)旋轉(zhuǎn)中心。運(yùn)用本文方法處理文獻(xiàn)[1]中的數(shù)據(jù)，能夠更精確地定位旋轉(zhuǎn)中心。

1 預(yù)備知識

1.1 投影正弦圖

在CT系統(tǒng)中，探測器-X射線源組合圍繞模板旋轉(zhuǎn)，一般為180°掃描，探測器接收到的信息即X射線經(jīng)過模板衰減后的強(qiáng)度。衰減后的強(qiáng)度受厚度、空氣、介質(zhì)的影響，在不同的投影視角，不同位置的探測器采集到投影數(shù)據(jù)，即旋轉(zhuǎn)角度θ與探測器t之間的關(guān)系為一條投影正弦曲線，即投影正弦圖。投影正弦圖中蘊(yùn)含著非常多的信息，其對稱性、正弦性、參數(shù)信息等都是不可忽略的重要內(nèi)容。

1.2 點(diǎn)擴(kuò)展函數(shù)與邊緣響應(yīng)函數(shù)

點(diǎn)擴(kuò)展函數(shù)（pointspreadfunction，PSF）在CT系統(tǒng)重建圖中表示空間中任意一點(diǎn)經(jīng)過CT系統(tǒng)投影、重建之后在圖像中的灰度強(qiáng)度分布[4-5]。邊緣響應(yīng)函數(shù)（edgeresponsefunction，ERF）是沿任意邊緣法線方向的灰度分布曲線。在重建圖像中，根據(jù)邊緣的任意一點(diǎn)法線方向灰度分布來繪制曲線，取任意方向上不同半徑的灰度值作曲線得邊緣響應(yīng)函數(shù)。物體的任意質(zhì)點(diǎn)經(jīng)過CT系統(tǒng)檢測和圖像重建之后，重建的圖像并非一個質(zhì)點(diǎn)，而是一個區(qū)域。獲取實(shí)際點(diǎn)擴(kuò)展函數(shù)的方法是，先獲取圖像的邊緣響應(yīng)函數(shù)，然后根據(jù)邊緣響應(yīng)函數(shù)計算圖像的點(diǎn)擴(kuò)展函數(shù)。

當(dāng)旋轉(zhuǎn)中心偏移量過大時，重建的圖像模糊，由單邊緣過渡到雙邊緣，邊緣出現(xiàn)一層穩(wěn)定灰度中間帶，點(diǎn)擴(kuò)展函數(shù)會出現(xiàn)雙峰，無法判斷極值與偏移量之間的關(guān)系。

1.3 Canny邊緣提取

在確定旋轉(zhuǎn)中心的過程中，需要通過Canny邊緣提取對投影圖、重建圖進(jìn)行處理。Canny邊緣檢測算法[6]首先用二維高斯函數(shù)的一階導(dǎo)數(shù)對圖像進(jìn)行平滑。

用分解的方法加快速度，把▽G的2個濾波卷積模板分解為2個一維行列濾波器，即：

式（2）中：k為常數(shù)；σ為濾波器參數(shù)，它控制著平滑程度；h1（t）＝ktexp（-t2/2σ2），h2（t）＝exp（-t2/2σ2）.

對于σ小的濾波器，雖然定位精度高，但信噪比低；而σ大的情況則相反。經(jīng)測試，本文高斯濾波器參數(shù)σ取0.27.

2 CT系統(tǒng)旋轉(zhuǎn)中心定位算法

2.1 確定初始旋轉(zhuǎn)中心

運(yùn)用PSF僅能精確在1個像素內(nèi)的偏移量[2]，對于偏移較遠(yuǎn)的旋轉(zhuǎn)中心，需要進(jìn)行初始定位。不同模板投影正弦圖有不同特征，本文僅考慮在橢圓模板情形下運(yùn)用投影正弦圖確定初始旋轉(zhuǎn)中心的問題。

對于N個探測器，從M個角度分別掃描得到投影圖，對任一像素位置 E（i，j）（i∈[0，M]，j∈[0，N]）做如下處理：①對于投影圖 E（i，j）（i∈[0，M]，j∈[0，N]），通過 Canny 邊緣提取，得到邊緣結(jié)果 T（i，j）（i∈[0，M]，j∈[0，N]）。②當(dāng)旋轉(zhuǎn)中心不發(fā)生偏移時，探測器中心與旋轉(zhuǎn)中心對齊，所以，可通過旋轉(zhuǎn)中心所在探測器位置來校正旋轉(zhuǎn)中心。通過T（i，j）計算當(dāng)前視角v下旋轉(zhuǎn)中心偏移量△Pv，單位為像素。

當(dāng)前視角v下旋轉(zhuǎn)中心偏移量△Pv，即為在該視角下偏移量的投影，如圖1所示，其中A為待定旋轉(zhuǎn)中心，B為精確的CT系統(tǒng)旋轉(zhuǎn)中心。

圖1 旋轉(zhuǎn)中心偏移量△Pv說明

橢圓模板的優(yōu)勢在于△Px、△Py對應(yīng)的投影數(shù)據(jù)容易判斷，探測器接收的數(shù)據(jù)最多，該列投影最長時對應(yīng)偏移量△Py；探測器接收的數(shù)據(jù)最少，該列投影最短對應(yīng)偏移量△Px，得到初始旋轉(zhuǎn)中心A＝（△Px，△Py）。

2.2 偏移旋轉(zhuǎn)中心計算及重建圖像

通過初始旋轉(zhuǎn)中心位置A或投影中心校正值d0校正旋轉(zhuǎn)中心。設(shè)校正后旋轉(zhuǎn)中心A1為p＝0的位置，PSF法僅能精確在1像素內(nèi)的偏移量，所以，從p＝1開始依次偏移0.1像素到p＝1和p＝-1的位置，并依次運(yùn)用直接反投影算法重建各個偏移量對應(yīng)的CT圖像。

2.3 基于Canny邊緣提取確定REF和PSF

傳統(tǒng)的Canny算法在2×2鄰域內(nèi)求有限差分來計算梯度幅值，該方法對噪聲比較敏感。本文采用在3×3鄰域內(nèi)計算梯度幅值[7]，考慮了像素對角線方向，增加了計算像素偏導(dǎo)數(shù)的方向，改進(jìn)了傳統(tǒng)梯度幅值計算方法，使得邊緣定位更準(zhǔn)確。建立邊緣響應(yīng)函數(shù)的具體步驟是：①在Matlab軟件中讀取重建的圖像，并通過rgb2gray命令得到灰度矩陣R＝（rij）n×n；②通過Canny邊緣提取結(jié)果確定出橢圓模板中心M（a，b），其半徑近似為r；③在模板中，以M（a，b）為起點(diǎn)，在水平正方向上分別取長度為r1和r2的2點(diǎn)，使得區(qū)間[r1，r2]跨越邊緣（本文取區(qū)間長度為30）。于是，邊緣響應(yīng)函數(shù)為x（j）＝raj，j∈[r1，r2].

在重建圖像中，對比度強(qiáng)、邊緣清晰的圖像邊緣響應(yīng)函數(shù)曲線的過渡帶短；對比度低、模糊的圖像，其邊緣響應(yīng)函數(shù)曲線的過渡帶相對比較長，可以跨越幾個像素甚至幾十個像素。邊緣響應(yīng)函數(shù)的一階導(dǎo)數(shù)為點(diǎn)擴(kuò)展函數(shù)y（j）。在工程實(shí)際中，點(diǎn)擴(kuò)展函數(shù)y（j）可通過差分格式計算，即y（j）＝Dx（j）＝x（j+1）-x（j），j∈[r1，r2].

2.4 CT系統(tǒng)旋轉(zhuǎn)中心的確定

點(diǎn)擴(kuò)展函數(shù)具有極大值，極值與偏移量間可以通過高斯函數(shù)擬合。當(dāng)旋轉(zhuǎn)中心不發(fā)生偏移，邊緣響應(yīng)函數(shù)陡峭變化，點(diǎn)擴(kuò)展函數(shù)極值最大。通過擬合點(diǎn)擴(kuò)展函數(shù)極值與偏移量，其高斯曲線峰值對應(yīng)的偏移量即為旋轉(zhuǎn)中心偏移量[2]。運(yùn)用PSF精確定位旋轉(zhuǎn)中心的步驟如下：①在不同偏移量下重建圖像，得到邊緣響應(yīng)函數(shù)和點(diǎn)擴(kuò)展函數(shù)；②求出不同偏移量下的點(diǎn)擴(kuò)展函數(shù)極值；③擬合點(diǎn)擴(kuò)展函數(shù)極值與偏移量之間的關(guān)系曲線G；④求出曲線G的峰值所對應(yīng)的偏移量，即旋轉(zhuǎn)中心偏移量。

由步驟①②③④分別得出x軸、y軸上偏移量△Px，△Py.如果以待確定旋轉(zhuǎn)中心的位置為原點(diǎn)，對應(yīng)的以重建圖像中心為原點(diǎn)建立直角坐標(biāo)系，則得到旋轉(zhuǎn)中心坐標(biāo)為（△Px，△Py）。

3 橢圓模板旋轉(zhuǎn)中心的確定

下面以2017年“高教社杯”全國大學(xué)生數(shù)學(xué)建模競賽A題給出的模板與相關(guān)數(shù)據(jù)為基礎(chǔ)確定其旋轉(zhuǎn)中心。

3.1 運(yùn)用投影正弦圖確定初始旋轉(zhuǎn)中心

賽題模板如圖2所示，根據(jù)接收信息畫出其投影正弦圖，如圖3所示。

先通過投影圖邊緣法確定初始旋轉(zhuǎn)中心，確定△Py＝22.5，△Px＝-33，得到初始旋轉(zhuǎn)中心A＝（-33，22.5）（單位為像素，下同）。依據(jù)初始旋轉(zhuǎn)中心A進(jìn)行校正，校正后的旋轉(zhuǎn)中心為A1，得到投影圖、重建圖如圖4、圖5所示。

3.2 偏移旋轉(zhuǎn)中心及重建圖像

水平方向為x方向，從校正后旋轉(zhuǎn)中心A1開始（設(shè)此時偏移量p＝0），分別依次偏移0.1像素到p＝1和p＝-1的位置，運(yùn)用直接反投影算法重建各個偏移量對應(yīng)的CT圖像。圖6給出了x方向上部分偏移量時的重建圖像。豎直方向（y方向）同理。

3.3 確定邊緣響應(yīng)函數(shù)

通過x（j）＝raj，j∈[r1，r2]得出不同偏移量下ERF函數(shù)。圖7給出了x方向、y方向上部分ERF曲線，從中可以看出其曲線有差異。旋轉(zhuǎn)中心偏移后越靠近系統(tǒng)精準(zhǔn)的旋轉(zhuǎn)中心，重建圖像越清晰，其ERF函數(shù)曲線陡峭下降，如圖7（a）中p＝-0.5，圖7（b）中p＝0.1所示；反之，重建圖像越模糊，其ERF函數(shù)曲線下降越平滑。

圖2 賽題模板示意圖

圖3 賽題模板投影正弦圖

圖4 校正后的投影圖

圖5 校正后的重建圖像

圖6 x方向不同偏移量時重建圖像

圖7 各方向上不同偏移量對應(yīng)的ERF函數(shù)曲線

3.4 確定點(diǎn)擴(kuò)展函數(shù)

通過 y（j）＝Dx（j）＝x（j+1）-x（j），j∈[r1，r2]得出不同偏移量下PSF函數(shù)。圖8給出x方向、y方向部分PSF曲線，從中可看出曲線有明顯差異。在圖8（a）中，從p＝-0.7到p＝-0.3，圖8（b）中從p＝-0.1到p＝0.3，對應(yīng)PSF極值先增大再減小，越接近系統(tǒng)精準(zhǔn)的旋轉(zhuǎn)中心，其PSF函數(shù)極值越大。

圖8 各方向上不同偏移量對應(yīng)的PSF曲線

3.5 根據(jù)旋轉(zhuǎn)中心偏移量和PSF極值確定旋轉(zhuǎn)中心

通過高斯函數(shù)擬合x方向、y方向上的偏移量與PSF函數(shù)極值的關(guān)系，得到圖9.圖9（a）峰值對應(yīng)x軸上偏移量為-0.527 3，圖9（b）峰值對應(yīng)0.103 6，旋轉(zhuǎn)中心A1＝（-0.527 3，0.103 6）；精確旋轉(zhuǎn)中心為A+A1＝（-33.527 3，22.603 6），單位為像素。

圖9 偏移量與PSF函數(shù)極值擬合

由PSF法精確后校正旋轉(zhuǎn)中心，再掃描模板并重建圖像。圖10由于旋轉(zhuǎn)中心偏移比較遠(yuǎn)，未校正重建圖像十分模糊，邊緣偽影、陰影嚴(yán)重；圖11初始校正后陰影減小，但邊緣依然不夠光滑；圖12通過PSF精確定位并校正后，重建圖像明顯清晰，邊緣光滑。

4 結(jié)束語

安裝CT系統(tǒng)會導(dǎo)致旋轉(zhuǎn)中心偏移，其誤差會使重建圖像出現(xiàn)偽影，降低成像質(zhì)量。系統(tǒng)旋轉(zhuǎn)中心位置定位越準(zhǔn)確，重建圖像質(zhì)量越好。本文先運(yùn)用投影正弦圖法進(jìn)行初始定位，再運(yùn)用PSF法進(jìn)一步精確定位了CT系統(tǒng)的旋轉(zhuǎn)中心。在未來的研究中，將用更具有特征的模板來標(biāo)定旋轉(zhuǎn)中心，例如均質(zhì)圓盤模板，其可以降低噪聲影響，便于分析和觀察。

圖12 PSF法校正

圖10 未校正

圖11 初始校正

［1］全國大學(xué)生數(shù)學(xué)建模競賽組委會.2017年競賽A題［EB/OL］.［引用日期不詳］.http://www.mcm.edu.cn.

［2］劉明進(jìn).工業(yè)CT系統(tǒng)旋轉(zhuǎn)中心定位方法研究［D］.重慶：重慶大學(xué)，2014.

［3］張蔚，羅守華，陳功.Micro-CT系統(tǒng)中對投影圖像旋轉(zhuǎn)中心的校正［J］.醫(yī)療衛(wèi)生裝備，2009，30（3）：4-6.

［4］孫怡，朱佩平，于健，等.X射線衍射增強(qiáng)成像中吸收、折射以及散射襯度的計算層析［J］.光學(xué)學(xué)報，2007，27（4）：749-754.

［5］張曉瑞.基于Radon變換的濾波反投影重建算法研究［J］.電腦知識與技術(shù)，2016，12（27）：259-261.

［6］王植，賀賽先.一種基于Canny理論的自適應(yīng)邊緣檢測方法［J］.中國圖象圖形學(xué)報，2004，9（8）：957-962.

［7］文章，張欣，周昌順，等.一種基于Canny的邊緣檢測改進(jìn)算法［J］.通信技術(shù)，2017，50（10）：2236-2240.

［8］孟凡勇，王維，李靜海.利用正弦圖自動校正CT投影中心［G］//全國射線數(shù)字成像與CT新技術(shù)研討會論文集.綿陽：中國體視學(xué)學(xué)會，2008：265-271.