基于數學建模下的CT圖形成像及系統參數標定

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

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

0 前言

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

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

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

1 CT系統工作與參數標定原理

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

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

圖1 CT系統示意圖

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

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

圖3 介質接收信息色階圖

幾何圖形位置情況如下：

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

圖4

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

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

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

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

圖5

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

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

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

圖6

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

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

圖7

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

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

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

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

圖8

坐標的參數方程可表示為

設正方形托盤中心為原點。

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

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

2 CT圖像模型的重建

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

■2.1 Radon變換

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

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

Radon 證明了以下定理：

圖9 Radon變換參數示意圖

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

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

clc,clear

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

R=radon(sheet1);

figure;

imshow(R,[]);

其中附件一里的數據是圖10的像素。

圖10 附件一讀出圖

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

圖11 radon變換后的圖

圖12 附件二讀出圖

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

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

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

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

第一步，設計合適的濾波器h。

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

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

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

圖13

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

clear all

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

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

imgsheet2cg=imresize(imgsheet2,256/362);

figure;imshow(imgsheet2cg,[]);

附件5 的數據是圖14的像素。

這段代碼的含義是實現radon反變換，所得圖像結果如圖15所示。

圖14

圖15 radon反變換的濾波圖

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

3 案例解析

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

圖16 10個位置示意圖

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

圖17 附件六濾波圖

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

圖18 附件六中的介質形狀圖

圖中，經過radon逆變換后，獲得圖像的吸收率矩陣ρ，觀察其形狀，我們可以看到圖由一些不規則的網格所構成，其吸收率最大值為以圖像矩陣的n/2進行柵格化處理，將其吸收率矩陣由離散值轉化為連續值，圖中白框為模板，紅圈為射線旋轉中心，藍點為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

再經過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

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

4 結論

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