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

巢麗媛，鄒娜，鐘仁峰，劉財(cái)輝*

?

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

巢麗媛，鄒娜，鐘仁峰，劉財(cái)輝*

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

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

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

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

引言

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

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

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

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

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

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

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

模型假設(shè)

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

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

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

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

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

符號約定與說明

（5）I1入射強(qiáng)度

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

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

（8）μ2附件二的吸收

2 問題一的分析與求解

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

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

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

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

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

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

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

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

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

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

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

2.2 問題一的求解

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 隨機(jī)選擇的兩條擬合圖像

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

（見附件擬合.xlsx）

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

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

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

表3 180 個(gè)角度值

(見附件 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 的線積分為：

此線積分表達(dá)式即為 Radon 變換。

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

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

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

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

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

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

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

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

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

3.2 模型二的建立

3.2.1 建立幾何圖形

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

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

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

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

3.2.2 獲取吸收率矩陣

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

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

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

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

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

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

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

圖8 附件 3 的物體信息

3.3.2 附件3物體的吸收矩陣

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

吸收率矩陣見附件 problem2.xls

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

十個(gè)點(diǎn)的坐標(biāo)以及對應(yīng)的吸收率：

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

4 問題三的分析與求解

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

圖11 附件 5 的物體信息

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

吸收率矩陣見附件 problem3.xls

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

十個(gè)點(diǎn)的坐標(biāo)及其對應(yīng)的吸收率：

表5 問題三對應(yīng)的十個(gè)點(diǎn)的吸收率

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

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

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

[4] 梁力, 尹東斐, 王川. 高精度攝像機(jī)標(biāo)定模板的設(shè)計(jì)及識(shí)別算法[J/OL]. 西安交通大學(xué)學(xué)報(bào), 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é)報(bào), 2015, 36(02): 71-74.

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

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

[12] 廖永忠, 蔡自興, 何湘華. 一種基于局部Radon變換運(yùn)動(dòng)模糊圖像參數(shù)估計(jì)算法[J]. 小型微型計(jì)算機(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ù)標(biāo)定及成像建模[J]. 數(shù)碼設(shè)計(jì), 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ù)計(jì)學(xué)院15 軟本2 班本科生；

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