CT系統(tǒng)參數(shù)標定及成像建模

巢麗媛，鄒娜，鐘仁峰，劉財輝*

?

CT系統(tǒng)參數(shù)標定及成像建模

巢麗媛，鄒娜，鐘仁峰，劉財輝*

（贛南師范大學(xué)數(shù)學(xué)與計算機科學(xué)學(xué)院，江西贛州，341000）

本文主要基于 Radon 正反變換，建模研究 CT 系統(tǒng)參數(shù)標定以及成像原理。針對問題1，以正方形模板的中心點為直角坐標系的原點，水平方向為 X 軸，豎直方向為 Y 軸，運用幾何方法、數(shù)據(jù)分析、圖像擬合等方法，解出CT 系統(tǒng)探測器單元之間的距離為0.2768mm，其旋轉(zhuǎn)中心的坐標為（-9.2734，5.9516）。CT系統(tǒng)的初始角度為33.0363°，終止角度為208.4130°，總共轉(zhuǎn)過 175.3767°。針對問題二和三，以問題一中得到的角度信息作為 Radon 反變換的參數(shù)，繪制附件 2 的 CT 系統(tǒng)成像圖。主要使用三次樣條插值與像素坐標轉(zhuǎn)化直角坐標等方法，解出成像圖的吸收率和幾何圖的吸收率之間的關(guān)系為：

μ1=2μ2，轉(zhuǎn)換得出幾何圖中各點吸收率值。

傅里葉函數(shù)；CFTOOL 工具箱；Radon 正反變換；三次樣條插值

引言

CT(Computed Tomography)可以在不破壞樣品的情況下，利用樣品對射線能量的吸收特性對生物組織和工程材料的樣品進行斷層成像，由此獲取樣品內(nèi)部的結(jié)構(gòu)信息。一種典型的二維 CT 系統(tǒng)如圖 1 所示，平行入射的 X 射線垂直于 CT 系統(tǒng)器平面，每個 CT 系統(tǒng)器單元看成一個接收點，且等距排列。X 射線的發(fā)射器和 CT 系統(tǒng)器相對位置固定不變，整個發(fā)射-接收系統(tǒng)繞某固定的旋轉(zhuǎn)中心逆時針旋轉(zhuǎn) 180 次。對每一個 X 射線方向，在具有 512 個等距單元的 CT 系統(tǒng)器上測量經(jīng)位置固定不動的二維待檢測介質(zhì)吸收衰減后的射線能量，并經(jīng)過增益等處理后得到 180 組接收信息。

CT 系統(tǒng)安裝時往往存在誤差，從而影響成像質(zhì)量，因此需要對安裝好的 CT 系統(tǒng)進行參數(shù)標定，即借助于已知結(jié)構(gòu)的樣品（稱為模板）標定 CT 系統(tǒng)的參數(shù)，并據(jù)此對未知結(jié)構(gòu)的樣品進行成像。

請建立相應(yīng)的數(shù)學(xué)模型和算法，解決以下問題：

（1）在正方形托盤上放置兩個均勻固體介質(zhì)組成的標定模板，模板的幾何信息如圖 2 所示，相應(yīng)的數(shù)據(jù)文件見附件 1，模板每一點的數(shù)值反映了該點的吸收強度，這里稱為“吸收率”。對應(yīng)于該模板的接收信息見附件 2。請根據(jù)這一模板及其接收信息，確定 CT 系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置、CT 系統(tǒng)器單元之間的距離以及該 CT 系統(tǒng)使用的 X 射線的 180 個方向。

（2）附件 3 是利用上述 CT 系統(tǒng)得到的某未知介質(zhì)的接收信息。利用(1)中得到的標定參數(shù)，確定該未知介質(zhì)在正方形托盤中的位置、幾何形狀和吸收率等信息。另外，請具體給出圖 3 所給的 10 個位置處的吸收率，相應(yīng)的數(shù)據(jù)文件見附件 4。

（3）附件 5 是利用上述 CT 系統(tǒng)得到的另一個未知介質(zhì)的接收信息。利用(1) 中得到的標定參數(shù)，給出該未知介質(zhì)的相關(guān)信息。另外，請具體給出圖 3 所給的10 個位置處的吸收率。

1 模型假設(shè)與符號約定

模型假設(shè)

（1）假設(shè) X 射線能夠全部覆蓋被測物體；

（2）將一條 X 射線抽象為一個質(zhì)點；

（3）假設(shè)每條 X 射線的入射強度不變。

（4）假設(shè)每個像素塊的吸收率是相同的

（5）假設(shè)所給束足夠可靠。

符號約定與說明

（5）I1入射強度

（6）μ物質(zhì)得衰減系數(shù)

（7）μ1附件一的吸收率

（8）μ2附件二的吸收

2 問題一的分析與求解

2.1 問題的分析及數(shù)據(jù)處理

問題一需要我們通過附件 1 和附件 2 以及模板的圖片等信息來標定一個 CT 系統(tǒng)的參數(shù)（包括旋轉(zhuǎn)中心，CT 系統(tǒng)單元之間的距離，以及系統(tǒng)旋轉(zhuǎn)的方向）。

分析附件二的圖像和系統(tǒng)旋轉(zhuǎn)的情況，可以得到系統(tǒng)旋轉(zhuǎn)時兩個特殊的方向

（水平方向和豎直方向），這兩個特殊方向的信息加上兩個物體的圖形信息，可以得到系統(tǒng)的旋轉(zhuǎn)中心和 CT 系統(tǒng)單元之間的距離。

兩種情況下，X 射線通過模板之后，所吸收的能量最大（如下圖所示）：

圖1 其中一條 X 射線穿過橢圓長軸

圖2 其中一條 X 射線穿過橢圓短軸

圖 1 是所有射線與水平方向成 90°的情況，此時，存在一條射線恰好通過 橢圓的長軸，在這個方向下，這條射線的 CT 系統(tǒng)單元返回的能量信息是所有射線中最大的。在附件二中篩選出這條射線的 CT 系統(tǒng)單元標號為 235，方向是第61 次旋轉(zhuǎn)所得方向。

圖 2 是所有射線與數(shù)值方向成 90°的情況，此時，存在一條射線恰好通過橢圓的短軸，在這個方向下，這條射線的 CT 系統(tǒng)單元返回的能量信息是所有射線中最大的，在附件二中篩選出這條射線的 CT 系統(tǒng)單元標號為 223，方向是第151 次旋轉(zhuǎn)所得方向。

表1 從附件 2 中篩選得到的 CT 系統(tǒng)單元與方向信息

（列標為方向，行標為 CT 系統(tǒng)單元標號）

2.2 問題一的求解

2.2.1 CT 系統(tǒng)單元之間的距離

在圖 3 的情況下，存在通過橢圓長軸的一簇射線束，這些射線束即為具有能量返回信息的射線束，可以在附件 2 中統(tǒng)計得到這一簇射線束中 X 射線的數(shù)量為289。X 射線的數(shù)量即為 CT 系統(tǒng)單元的數(shù)量 289。

通過題中給出的模板示意圖，已知條件容易得到橢圓長軸的長度為 80mm。80mm對應(yīng) 289 個 CT 系統(tǒng)單元的 288 段間距，而每段間距又是相等的，故使用d表示 CT 系統(tǒng)單元之間的距離，a 表示橢圓長軸，用如下公式計算。

求得 CT 系統(tǒng)單元之間的距離d = 0.2768 （單位：mm）。

2.2.2 旋轉(zhuǎn)中心在正方形托盤中的位置

在整個 CT 系統(tǒng)中，旋轉(zhuǎn)中心的位置一定處于 CT 系統(tǒng)單元 256 與 CT 系統(tǒng)單元 257 所發(fā)射出的兩條 X 射線的中間。

為了確定旋轉(zhuǎn)中心在正方形托盤中的位置，我們需要將 CT 系統(tǒng)與正方形托盤上的模板聯(lián)系起來看，當 CT 系統(tǒng)旋轉(zhuǎn)到 3.1 中圖1 和圖2 所示的兩個位置時，分別有兩個 CT 系統(tǒng)單元（235、223）發(fā)出的射線恰好經(jīng)過橢圓或圓的中心。CT 系統(tǒng)單元 256、257 與這兩個 CT 系統(tǒng)單元（235、223）的相對位置可以確定，在圖像中心建立笛卡爾直角坐標系，就可以得到旋轉(zhuǎn)中心的坐標。

計算吸收量最大的 CT 系統(tǒng)單元與 CT 系統(tǒng)單元 256 和 CT 系統(tǒng)單元 257 之間的間隔，間隔乘以相鄰 CT 系統(tǒng)單元之間的距離，即可得到旋轉(zhuǎn)中心距離 x 軸與y 軸的距離（假設(shè) 235 號射線與 223 號射線是 x 軸和 y 軸，橢圓的中心是原點 O）。

圖3 旋轉(zhuǎn)中心在圖中的位置

旋轉(zhuǎn)中心的坐標（-9.2734，5.9516）（單位：mm）

2.2.3 旋轉(zhuǎn)180次的方向

二維圖像重建算法常用傅里葉變換。設(shè) f (x, y) 表示一個二維圖像，其傅里葉變換為：

這個二維圖像在與 x 軸夾? 角的射線 s 上的投影的傅里葉變換，恰好等于這個二維圖像 f (x, y) 的二維傅里葉變換[3]。

CT 系統(tǒng)旋轉(zhuǎn)時，二維圖像 x 軸夾θ角的射線 s 上的投影在變化，也就是傅里葉函數(shù)在變化，這些變化恰好能體現(xiàn)這 180 次旋轉(zhuǎn)方向的變化，我們認為這種角度的變化滿足傅里葉函數(shù)。在附件 2 中，隨機選取幾條 X 射線所吸收的能量值隨 180 個方向的變化，使用 MATLAB 中 CFTOOL 工具箱擬合圖像，如圖4所示。

圖4 隨機選擇的兩條擬合圖像

表2 兩條曲線擬合的結(jié)果

（見附件擬合.xlsx）

上圖和上表容易看出，擬合的精度很高。

將擬合所得的結(jié)果處理之后，得到角度變化所滿足的傅里葉函數(shù)：

將得到的曲線的函數(shù)值（即對應(yīng)附件 2 中的信息）帶入上述函數(shù)，解得對應(yīng)f(x)下x的值，即 180 個變化的θ值。

表3 180 個角度值

(見附件 theta.xlsx)

上述數(shù)據(jù)中，初始角度0 ?33.0363，終止角度t? 208.4130?，轉(zhuǎn)過的角度175.3767.

3 問題二的模型建立與求解

3.1 問題二的數(shù)據(jù)分析

3.1.1 Radon 變換與反變換

定義函數(shù) f (x, y) 在平面上沿直線 L 的線積分為：

此線積分表達式即為 Radon 變換。

數(shù)學(xué)家 Radon 給出了上述公式 4.1 的逆變換表達式：

Q 是坐標軸上任意一點，q 是直線 L 到點 Q 的距離， FQ (q) 是關(guān)于 L 的線積分 Pf (L) 對所有 q 的平均值，公式 4.2 即為 Radon 反變換。

根據(jù) Beer-Lamber 定律，X 射線的強度衰減滿足：

當X 射線穿過不同衰減系數(shù)的材料組成的非均勻物體時，公式4.3 應(yīng)該寫成：

其中?(,)是 ?(,) 沿 L 的線積分，化簡公式 4.4，會得到一個類似于Radon 變換的形式：

3.1.2 利用標定參數(shù)得到幾何圖

問題二需要我們求解附件 3 中未知介質(zhì)在正方形托盤中的位置、幾何形狀和吸收率等信息。首先我們需要在 CT 系統(tǒng)成像圖中得到原幾何圖的信息。通過問題 1 中標定得到的參數(shù)——180 個旋轉(zhuǎn)角度以及附件 2 中的能量信息，使用 iradon 函數(shù)繪制出模板的 CT 系統(tǒng)成像圖像。這一點容易用 MATLAB 實現(xiàn)。

通過問題 1 中標定得到的參數(shù)——旋轉(zhuǎn)中心和 CT 系統(tǒng)單元之間的距離，在CT 系統(tǒng)成像圖中選取出原幾何圖的信息。

根據(jù)計算，得到 CT 系統(tǒng)成像圖的對角線寬是 512 個 CT 系統(tǒng)單元的總寬度，即 CT 系統(tǒng)單元數(shù)與 CT 系統(tǒng)單元之間距離的乘積。而且， CT 系統(tǒng)成像圖的中心是 CT 系統(tǒng)的旋轉(zhuǎn)中心。如此，可以標定 CT 系統(tǒng)成像圖的實際坐標（單位： mm），通過實際坐標與原幾何圖標定坐標之間的關(guān)系確定原幾何圖的中心和幾何圖的邊緣。

3.2 模型二的建立

3.2.1 建立幾何圖形

首先，將通過 radon 變換得到的像素為 362 的圖像擴展為像素 512 的圖像，此時圖像的實際長度為 141.4447mm。

在 CT 系統(tǒng)圖的中心選擇坐標為（256.5， 256.5），即為旋轉(zhuǎn)中心的坐標，對應(yīng)原幾何圖中的實際坐標為（-9.2734， 5.9516）（單位： mm）。

利用問題一中 CT 系統(tǒng)器水平方向與豎直方向的 X 射線 223 和 235 與中心線之間的關(guān)系，幾何圖中心點在 CT 系統(tǒng)成像圖上的坐標為（289， 277）。關(guān)于幾何圖的邊緣計算，有如下等式：

解得 x=y=361,即白色邊框的邊長。

3.2.2 獲取吸收率矩陣

圖6 幾何圖吸收率的三維圖像

圖7 CT系統(tǒng)圖吸收率的三維圖像通過

幾何圖吸收率三維圖像與CT系統(tǒng)圖吸收率三維圖像的吸收率之間的關(guān)系得到：

由于在 CT 系統(tǒng)成像圖中得到的幾何圖的像素圖大小是 361*361 的矩陣，而幾何圖的本身像素大小是 256*256 的矩陣，為了實現(xiàn)成像圖中幾何圖信息到實際幾何圖信息的轉(zhuǎn)換，我們對 CT 系統(tǒng)圖吸收率的三維圖像進行插值處理，得到更詳細的吸收率值。

使用插值后的成像圖，根據(jù)幾何圖和成像圖的邊長比，將幾何圖的距離映射到插值后的成像圖中，即求出 256*256 像素塊在插值后的成像圖中的坐標，根據(jù)幾何圖吸收率與成像圖的上述關(guān)系，從而求出幾何圖中每個像素塊的吸收率。

3.3 模型二在附件 3 的應(yīng)用

3.3.1 附件 3 物體在正方形托盤中的位置，以及幾何形狀

圖8 附件 3 的物體信息

3.3.2 附件3物體的吸收矩陣

圖 9 附件 3 的吸收率三維圖像

吸收率矩陣見附件 problem2.xls

圖10 通過吸收率矩陣重建的圖像

十個點的坐標以及對應(yīng)的吸收率：

表4 問題二對應(yīng)的十個點的吸收率值

4 問題三的分析與求解

與模型二中的處理方法類似，首先得到 CT 系統(tǒng)成像圖，在 CT 系統(tǒng)成像圖中標定幾何圖的位置，如下圖中白色方框框選的部分。

圖11 附件 5 的物體信息

圖12 附件 5 的吸收率的三維圖像

吸收率矩陣見附件 problem3.xls

圖13 通過吸收率矩陣重建的圖像

十個點的坐標及其對應(yīng)的吸收率：

表5 問題三對應(yīng)的十個點的吸收率

[1] 姜啟源, 謝金星, 葉俊著. 數(shù)學(xué)模型（第四版）[M]. 高等教育出版社, 2011.

[2] 阮秋琦著. 數(shù)字圖像處理學(xué)[M]. 電子工業(yè)出版社, 2013.

[3] 朱翚, 王富東. 利用 MATLAB 實現(xiàn)二維圖像傅立葉變換算法[J]. 計算機應(yīng)用與軟件, 2006, (12): 141-142.

[4] 梁力, 尹東斐, 王川. 高精度攝像機標定模板的設(shè)計及識別算法[J/OL]. 西安交通大學(xué)學(xué)報, 2011, 45(04): 82-85.

[5] 李翰威. 錐形束CT系統(tǒng)幾何偽影校正技術(shù)研究[D]. 南方醫(yī)科大學(xué), 2015.

[6] 陳健, 黃政仁, 劉學(xué)建, 等. CT系統(tǒng)放大倍數(shù)與極限空間.

[7] 姚明, 叢鵬, 劉錫明. 平板探測器CT系統(tǒng)調(diào)校及參數(shù)獲取方法[J]. 原子能科學(xué)技術(shù), 2013, 47(06): 1019-1022.

[8] 劉興龍, 孫宏, 李琛瑋, 等. CT能量成像技術(shù)原理和應(yīng)用[J]. 中國醫(yī)療設(shè)備, 2012, 27(09): 1-7+45.

[9] 劉文艷, 許友軍, 劉亞春. MATLAB在傅里葉變換課堂教學(xué)中的應(yīng)用[J]. 湘南學(xué)院學(xué)報, 2015, 36(02): 71-74.

[10] 康家方, 王紅星, 鐘佩琳, 等. 基于傅里葉變換的非正弦時域正交調(diào)制系統(tǒng)實現(xiàn)方法[J]. 上海交通大學(xué)學(xué)報, 2014, 48(10): 1415-1420.

[11] 畢清華. Radon變換奇性檢測與反演的研究[D]. 北京交通大學(xué), 2007.

[12] 廖永忠, 蔡自興, 何湘華. 一種基于局部Radon變換運動模糊圖像參數(shù)估計算法[J]. 小型微型計算機系統(tǒng), 2014, 35(01): 133-136.

[13] 王麗新. 雙能CT基物質(zhì)分解算法應(yīng)用研究[D]. 山東大學(xué), 2016.

[14] 桂葉晨. 基于CUDA的錐束CT重建與CT圖像可視化技術(shù)研究[D]. 南方醫(yī)科大學(xué), 2009.

[15] CT成像: 基本原理、偽影與誤區(qū)[J]. 中國介入影像與治療學(xué), 2015, 12(08): 492.

Parameter Calibration and Mmodeling of CT System

CHAO Liyuan, ZOU Na, ZHONG Renfeng, LIU Caihui*

(School of mathematics and computer science, Gannan Normal University, Jiangxi Ganzhou, 341000, China)

Fourier Transform; CFTOOL; Radon Transform; Cubic Spline Interpolation

10.19551/j.cnki.issn1672-9129.2018.01.012

TP39

A

1672-9129(2018)01-0027-05

巢麗媛, 鄒娜, 鐘仁峰, 等. CT系統(tǒng)參數(shù)標定及成像建模[J]. 數(shù)碼設(shè)計, 2018, 7(1): 27-31.

CHAO Liyuan, ZOU Na, ZHONG Renfeng, et al. Parameter Calibration and Mmodeling of CT System[J]. Peak Data Science, 2018, 7(1): 27-31.

2017-10-21；

2017-12-13。

國家自然科學(xué)基金(61305052)，贛南師范大學(xué)校級教改課題(150656)。

巣麗媛（1997-），女，贛南師范大學(xué)數(shù)計學(xué)院15 軟本2 班本科生；

劉財輝（1979-），男，副教授，博士，碩士生導(dǎo)師，研究方向：數(shù)據(jù)智能處理與應(yīng)用、機器學(xué)習(xí)等。E-mail:839783509@qq.com