基于數(shù)學(xué)建模下的CT圖形成像及系統(tǒng)參數(shù)標定

王美晨，李嘉萱，席典兵，孫新宇

（1.南昌大學(xué)信息工程學(xué)院，江西南昌，330031；2.南昌大學(xué)資源環(huán)境與化工學(xué)院，江西南昌，330031；3.南昌大學(xué)材料科學(xué)與工程學(xué)院，江西南昌，330031；4.南昌大學(xué)信息工程學(xué)院，江西南昌，330031）

0 前言

CT(Computed Tomography)，即電子計算機斷層掃描，具有圖像清晰、掃描時間少等特點，它通過利用精確準直的γ射線、X線束、超聲波等，與靈敏度極高的探測器一同圍繞人體的某一部位作一個接一個的斷面掃描[1]。 CT技術(shù)應(yīng)用廣泛，可用于多種疾病的檢查也可應(yīng)用于重要產(chǎn)品的檢測以及材料疲勞、汽車鑄件、裝配分析、焊縫三維成像等眾多領(lǐng)域[2]。

文獻[3]提出利用雙層字典學(xué)習(xí)方法對圖像進行重建，與傳統(tǒng)的K-SVD算法相比較雙層字典學(xué)習(xí)方法得出的圖像更接近于原始CT圖像，提升了圖像的質(zhì)量。文獻[4]介紹了對CT圖像后期處理技術(shù)的評價方法，通過對兩種技術(shù)后期臨床診斷效果的比較，得出16層螺旋CT圖像后處理技術(shù)比常規(guī)組效果顯著。文獻[5]通過一種指數(shù)形式的濾波函數(shù)(EBFF)對圖像進行重建，與傳統(tǒng)的濾波函數(shù)相比較，基于EBFF下重建的圖像質(zhì)量更優(yōu)，更接近于原始的CT圖像。文獻[6]介紹了CT圖像在集料棱角性的評估上的應(yīng)用。文獻[7]介紹了如何利用流處理模型，加速了錐-束CT圖像的重建和三維變形圖像的使用。

近些年來國內(nèi)外對CT圖像的研究大多是對CT圖像重建的算法進行優(yōu)化或者是CT圖像技術(shù)在相關(guān)領(lǐng)域的應(yīng)用研究。很少文獻對CT圖像的成像原理以及CT系統(tǒng)參數(shù)的標定進行深入詳細的分析。而CT圖像的成像原理以及系統(tǒng)參數(shù)的標定是一切有關(guān)CT技術(shù)研究的基石。故本文主要研究CT系統(tǒng)參數(shù)標定及成像原理，對CT系統(tǒng)的幾何成像進行了深入淺出的分析。

1 CT系統(tǒng)工作與參數(shù)標定原理

在深入介紹CT圖像重建之前，我們先利用一個簡單的幾何體去介紹CT系統(tǒng)工作原理。以2017年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽題目為例，CT可以利用樣品對射線能量的吸收特性對生物組織和工程材料的樣品進行斷層成像，并且在不破壞樣品的情況下獲取樣品內(nèi)部的結(jié)構(gòu)信息。圖1是一種二維平行束CT系統(tǒng)， X射線平行入射并且垂直于探測器所在的平面，每個探測器等距排列，并且探測器單元可看成一個接收點。

假設(shè)在正方形托盤上放置兩個由均勻固體介質(zhì)組成的標定模板，標定模板如圖2所示。當標定模板工作在CT系統(tǒng)中時。我們需根據(jù)這一模板和其接收信息，確定探測器單元之間的距離、CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的具體位置以及該CT系統(tǒng)使用的X射線的180個方向。X射線的發(fā)射器和探測器相對位置固定不變時，整個發(fā)射-接收系統(tǒng)繞某固定的旋轉(zhuǎn)中心逆時針旋轉(zhuǎn)180次。2017年全國大學(xué)生數(shù)學(xué)建模大賽A題附件二便是此時對每一個X射線方向，在具有512個等距單元的探測器上測量經(jīng)位置固定不動的二維待檢測介質(zhì)吸收衰減后的射線能量，并經(jīng)過增益等處理后得到的180組接收信息 。附件2中每一點的數(shù)值體現(xiàn)了該點的吸收強度，也可以被稱作為“吸收率”[8]。

圖1 CT系統(tǒng)示意圖

圖2 模板示意圖（單位：mm）

用excel對附件2進行色階處理，根據(jù)題意并結(jié)合色階處理圖分析，其實被色階處理過的附件2中的接受信息可以看作是橢圓與小圓的幾何投影。圖3部分之所以有深有淺是因為接收信息不同，本模型的介質(zhì)是均勻的，所以吸收信息能力是相同的，所以理解吸收數(shù)據(jù)信息的不同可以轉(zhuǎn)化為被X射線照射時的厚度不同，簡而言之，某塊地方接收數(shù)據(jù)信息大，則說明這塊被X射線照射的地方厚。可以得出圖中的帶狀部分為標定模板上的圓形固體介質(zhì)的接受信息，剩下的瓶狀部分為標定模板上的橢圓形固體介質(zhì)的接受信息。對以下色階圖進行分析我們可以獲得五列特殊方位的數(shù)據(jù)。

圖3 介質(zhì)接收信息色階圖

幾何圖形位置情況如下：

（1）當附件2中數(shù)據(jù)處于第61列時（X射線方向如圖4所示）

圖4

當角度轉(zhuǎn)為以上情況時，在附件2中小圓和橢圓色階必定重合，被橢圓形固體介質(zhì)所吸收的X射線能量所占的跨度最大，并且無論哪個角度小圓形固體介質(zhì)所吸收的X射線能量所占跨度都是一樣的，即可通過分析附件2表格中哪一列的非零數(shù)最多，判斷出瓶狀圖形左側(cè)最粗的位置所占的列，則該列數(shù)據(jù)所表示的X射線的方向為上圖所示。

由MATLAB求解可得，附件2中第 58、59、60、61、62、63、64、65列所含的非零元素最多，即可從第61列和第62列中選出一列作為最佳上述情況。由于第61列中最大的數(shù)據(jù)大于第62列中最大的數(shù)據(jù)，即選擇第61列為最佳列。

（2）當附件2中數(shù)據(jù)處于第151列時（X射線方向如圖5所示）

當角度轉(zhuǎn)為以上情況時，被橢圓形固體介質(zhì)所吸收的X射線能量所占的跨度最小，且無論哪個角度小圓形固體介質(zhì)所吸收的X射線能量所占跨度都是一樣的，即可通過分析附件2表格中哪一列的非零數(shù)個數(shù)最少，判斷出瓶狀圖形右側(cè)最細的位置所占的列，則該列數(shù)據(jù)所表示的X射線的方向為圖5所示。

圖5

由MATLAB求解可得，第150、152列所含非零元素最少。但是發(fā)現(xiàn)最大的數(shù)據(jù)出現(xiàn)在第151列中，所以在誤差允許范圍之內(nèi)，我們選擇第151列作為最佳列。

由圖4至圖5可明確發(fā)現(xiàn)整個發(fā)射-接收系統(tǒng)繞某固定點逆時針旋轉(zhuǎn)了90°，且獲得了90（151-61）組數(shù)據(jù)，由此可以得出旋轉(zhuǎn)90°可以獲得90組接收信息的結(jié)論，進一步推斷為獲得180組接收信息需要旋轉(zhuǎn)180°。對帶狀數(shù)據(jù)的進一步擬合可以發(fā)現(xiàn)該帶狀數(shù)據(jù)可近似看做一正弦函數(shù)，即可得整個發(fā)射-接收系統(tǒng)是繞某固定點勻速旋轉(zhuǎn)的，即每旋轉(zhuǎn)一度記錄一組數(shù)據(jù)。

（3）當附件2中數(shù)據(jù)處于第0列也就是初始位置時（X射線方向如圖6所示θ=-61°）

圖6

圖6表示的是初始位置，由之前的結(jié)論可得當出現(xiàn)圖4時，X射線已逆時針旋轉(zhuǎn)了61°。現(xiàn)如上圖所示建立極坐標，坐標軸為圖示中的水平向右的加粗箭頭，易得出初始位置為θ=-61°。以下所有圖均以此方向建立極坐標軸。

（4）當附件2中數(shù)據(jù)處于第145列時（X射線方向如圖7所示θ=84°）

圖7

當X射線角度轉(zhuǎn)為以上情況時，所對應(yīng)的列數(shù)為第145列，對應(yīng)附件2中正弦函數(shù)圖像的最高點所對應(yīng)的列，此時旋轉(zhuǎn)中心與小圓圓心的連線應(yīng)垂直于X射線。所以可推斷出旋轉(zhuǎn)中心在板的上半部分。

（5）當附件2中數(shù)據(jù)處于第55列時（X射線方向如圖8所示θ=-6°）

該列數(shù)值為圖7中的X射線以順時針方向往回轉(zhuǎn)90°所得，此時旋轉(zhuǎn)中心與小圓圓心的連線是平行于X射線方向的。

令旋轉(zhuǎn)中心到小圓圓心的距離為半徑R，可通過求得第145列正弦條狀數(shù)據(jù)中最大值對應(yīng)的行數(shù)為59，以及第55列正弦條狀數(shù)據(jù)中最大值對應(yīng)的行數(shù)為255，它們之間的行差數(shù)乘以單元格間的距離即為半徑R，

圖8

坐標的參數(shù)方程可表示為

設(shè)正方形托盤中心為原點。

由MATLAB編程，可解得旋轉(zhuǎn)點坐標為（-8.8140，5.6561）。

該CT系統(tǒng)使用X射線的180個方向可表示為θ=-61+n（n=0,1,2…179）。此外此模型還可標定CT系統(tǒng)的參數(shù)。CT系統(tǒng)在初次使用時需要進行參數(shù)標定，避免產(chǎn)生誤差影響成像的質(zhì)量。

2 CT圖像模型的重建

由以上的幾何模型的建立可以說是為我們理解CT圖像的重建打下了堅實的基礎(chǔ)，自此開始我們將深入講述CT圖像重建的原理。

■2.1 Radon變換

奧地利數(shù)學(xué)家Radon在1917年提出了投影圖像重建的基本數(shù)學(xué)理論，為CT技術(shù)建立了數(shù)學(xué)理論基礎(chǔ)。

該理論從數(shù)學(xué)上證明了某種物理參量(如一個切片衰減系數(shù)的分布)的二維分布函數(shù)，由該函數(shù)在其定義域內(nèi)所有線積分完全確定。如圖9所示，二維平面內(nèi)的一條直線L與X軸夾角為?，原點到L的垂線距離為s，直線上的點(x , y )可以用極坐標表示為( r ,θ)[2]。

Radon 證明了以下定理：

圖9 Radon變換參數(shù)示意圖

式(1)被稱為Radon變換，是指二維分布函數(shù)在一定角度下的線積分，即實際的射線投影 p。式(2)被稱為Radon反變換，Radon反變換對CT的重建有重要的理論意義，它指的是通過一定量的投影采樣角度下的投影數(shù)據(jù)p可以重建出物體的斷層圖像

這個原理看上去很復(fù)雜，但當我們用代碼實現(xiàn)的時候便簡單明了了許多，Matlab有自帶的函數(shù)可以實現(xiàn)Radon變換。函數(shù)代碼如下[9]：

clc,clear

sheet1=xlsread(‘附件.xls’,’附件1’);

R=radon(sheet1);

figure;

imshow(R,[]);

其中附件一里的數(shù)據(jù)是圖10的像素。

圖10 附件一讀出圖

這段代碼的含義實際是將圖10進行radon變換，變換結(jié)果如圖11所示。

圖11 radon變換后的圖

圖12 附件二讀出圖

這幅圖與上文幾何模型提及的附件二的圖驚人地相似(圖12)，不同點是因為Matlab默認radon變換是以幾何物體中心為坐標原點變換的，與我們附件二的坐標原點不相同。

■2.2 平行束濾波反投影重建算法及其計算機實現(xiàn)流程圖

濾波反投影(Filtered back-projection, FBP)重建算法是最普遍的CT重建算法之一。待重建的圖像為由傅立葉中心切片定理推導(dǎo)而來的濾波反投影重建公式為：

由上式可知，濾波反投影重建圖像的具體過程為：首先把線陣探測器上獲得的投影數(shù)據(jù)與濾波器函數(shù)進行卷積運算，得到各方向卷積濾波后的投影數(shù)據(jù)；然后再沿各方向進行反投影，即依照其原路徑把它們平均分配到每一矩陣單元上，進行疊加后得到每一矩陣單元的CT值；進行些許適當處理，最后被掃描物體的斷層圖像完成[1]。算法步驟為：

第一步，設(shè)計合適的濾波器h。

第二步，把濾波器h與恒定角度Φi情況下所測投影進行卷積濾波，得到濾波后的投影

第三步，對于每一個Φi，把濾波后的投影反投影于符合公式的射線上的所有各點(r ,θ)。

第四步，將步驟三中所有0≤Φi≤2π的反投影值累加，從而得出重建后圖像數(shù)據(jù)圖13[1]。

圖13

此原理簡單地說就是radon變換的逆推導(dǎo),與radon變換一樣，Matlab有自帶的函數(shù)可以實現(xiàn)此原理。函數(shù)代碼如下[9]：

clear all

sheet2=xlsread(‘附件.xls’,’附件5’);

imgsheet2=iradon(sheet2,[0:179],’nearest’,’Ha mming’);

imgsheet2cg=imresize(imgsheet2,256/362);

figure;imshow(imgsheet2cg,[]);

附件5 的數(shù)據(jù)是圖14的像素。

這段代碼的含義是實現(xiàn)radon反變換，所得圖像結(jié)果如圖15所示。

圖14

圖15 radon反變換的濾波圖

我們由圖10發(fā)現(xiàn)圖片上的橢圓是斜的，這是因為上文中圖6可知，初始時刻的角度θ=?61°。

3 案例解析

假設(shè)利用幾何模型的CT系統(tǒng)得到的另一個未知介質(zhì)的接收信息附件六。利用上文中得到的標定參數(shù)，給出該未知介質(zhì)的相關(guān)信息。另外，請具體給出圖16所給的10個位置處的吸收率。

圖16 10個位置示意圖

在轉(zhuǎn)化為濾波圖的過程中，在讀入的矩陣前后添加zeros（100,180），使后期圖像位置更為美觀與處理，如圖17所示。

圖17 附件六濾波圖

通過上文可知，我們得到標定初始角度為-61度，以此為初始角度，進行radon反變換，獲得圖像形狀及其在標定模板的位置。獲得圖像如圖18所示。

圖18 附件六中的介質(zhì)形狀圖

圖中，經(jīng)過radon逆變換后，獲得圖像的吸收率矩陣ρ，觀察其形狀，我們可以看到圖由一些不規(guī)則的網(wǎng)格所構(gòu)成，其吸收率最大值為以圖像矩陣的n/2進行柵格化處理，將其吸收率矩陣由離散值轉(zhuǎn)化為連續(xù)值，圖中白框為模板，紅圈為射線旋轉(zhuǎn)中心，藍點為10個待求點，其坐標以（50-xc,50+yc）進行變換，坐標位置如下：

坐標點 橫坐標 縱坐標 坐標點 橫坐標 縱坐標1 -48.9713 -37.356 6 -8.9713 20.144 2 -24.4713 -30.356 7 -2.9713 21.144 3 -15.4713 -22.356 8 6.5287 -18.356 4 -13.9713 20.144 9 20.5287 -37.356 5 -10.4713 0.144 10 39.5287 -11.856

再經(jīng)過matlab程序求解，可得出10個位置。

坐標點 吸收率 坐標點 吸收率1-0.0067 6 1.0478 2 0.0114 7 1.0548 3 0.0042 8 0.0104 4 1.0442 9 0.0104 5 1.0258 10 0.0357

與上圖中的圖像進行對比，吻合的很好，最后再將其導(dǎo)出至model.xls。

4 結(jié)論

本文結(jié)合了圖片、實例生動形象地介紹了CT圖像的重建原理，其中還涉及CT系統(tǒng)參數(shù)標定及成像問題，這需要先借助模板來標定CT系統(tǒng)參數(shù)，再將參數(shù)運用于確定未知介質(zhì)的位置、幾何形狀和吸收率等信息。案例結(jié)果表明，本文提出的方法有一定的可行性。