基于滤波反投影重建成像及CT系统参数标定

荣亚亚 郭子健 吉 督 赵旭东

(1.华北理工大学经济学院 河北 唐山 063210;2.华北理工大学理学院 河北 唐山 063000)(3.华北理工大学机械工程学院 河北 唐山 063210;4.华北理工大学建筑工程学院 河北 唐山 063210)

基于滤波反投影重建成像及CT系统参数标定

荣亚亚1郭子健2吉 督3赵旭东4

(1.华北理工大学经济学院 河北 唐山 063210;2.华北理工大学理学院 河北 唐山 063000)(3.华北理工大学机械工程学院 河北 唐山 063210;4.华北理工大学建筑工程学院 河北 唐山 063210)

CT即计算机断层成像技术，是对物体进行不同角度的扫描投影重建而获得物体截面信息的成像技术。本文利用X射线、探测器与标定模板几何信息对一维CT系统进行研究，通过滤波反投影重建算法对断层图像进行处理，实现了滤波反投影图像重建，给出了一种新型的标定模板。

R－L滤波器;Radon逆变换;反投影图像重建

针对问题一，首先，解决探测器单元的距离问题:以标定模板的椭圆中心为坐标原点，椭圆圆心和小圆圆心连线所在直线为X轴建立直角坐标系;对探测器接收信息数据矩阵进行分析，建立了小圆实际半径与小圆投影的探测器接收信息数据矩阵关系，求解得出探测器单元之间距离为0.2808mm，总长为143.4388mm;其次，解决X射线的180个方向问题:建立了基于椭圆方程和欧式距离公式等几何关系的转角计算模型，求解出转角中X射线的初始位置约为－58.6900°;最后，解决CT系统的旋转中心问题:当X射线分别平行于X轴和Y轴时，对探测器中心，两圆中心及坐标轴进行几何分析，确定旋转中心坐标为 (－9.5356， －6.3180)。

针对问题二，首先，对探测器的接收信息矩阵数据进行分析处理，选取各个方向的投影;其次，用R－L滤波器对数据进行滤波处理，建立Radon逆变化模型，进行反投影图像重建，对重建的图像分析，通过坐标系转换对重建图像进行调整修复;最后，确定出未知介质在正方形托盘中几何位置及几何形状，对重建图像矩阵进行像素压缩，压缩为256*256的图像矩阵，建立压缩后的矩阵与正方形托盘中实际坐标的映射关系，得出正方形托盘上任意位置的吸收率。针对题中所给10个位置，求解出吸收率分别为0.0012，0.5103，－0.0002，0.4912，0.4936，0.4932，0.4952，0.5585， －0.0027， －0.0021。

针对问题三，对未知介质的接收信息矩阵，建立基于问题二的滤波反投影图像重建模型，采取相同技术处理路线，确定此未知介质在正方形模板的几何位置，几何形状和吸收率信息。计算得出10个位置的吸收率分别为 0.3976，0.4943，0.7050， －0.0424，1.2388，0.5597， －0.0738，1.0336，0.2364，0.2505。

针对问题四，首先，对原始模板的参数标定进行精度分析，结果表明探测器单元间距离误差为0.0015、CT系统旋转中心在X，Y轴的坐标误差为0.4050，0.3077。创新模板采用菱形与圆形组合的模板与矩形模板，对探测器单元距离和旋转中心参数进行标定，提高了原始模板的参数标定值精度;其次，分别对参数标定的稳定性和系统稳定性进行优化分析，对原始模板存在参数标定时会有波动范围和吸收率的噪声点的出现;最后，针对噪声点给出解决方案:引入反馈调节，对噪声点数据进行识别处理，选取合适的模板材料。最终，模拟绘制出3种创新组合模板和对应标定模型。

一、问题分析

问题一的分析:确定以及该CT系统使用的X射线的180个方向。首先，由于小圆在探测器上的投影是不变的，根据小圆的直径与投影的关系计算得出探测器单元之间的距离。CT系统旋转中心在正方形托盘中的位置的确定转化成特殊位置求解，垂直于X轴和垂直于Y轴两个位置。再根据实际距离和探测器单元之间的距离的关系求出旋转中心的坐标。

问题二的分析:利用问题一得到的标定参数，解决确定未知介质在托盘中的位置，几何形状，吸收率和托盘上10个位置点的吸收率。使用滤波反投影图像重建的方法，采用的技术路线为:(1)取投影 (2)数据滤波处理 (3)坐标系变换 (4)运用Radon逆变换 (5)图像反投影重建。使用了Matlab等相关工具包进行编码绘制图像的几何形状，在坐标系中的几何位置，并计算出180次选择，各个位置处的吸收率。

问题三的分析:基于问题二，给出未知介质的几何形状，几何位置，吸收率，对问题三的数据进行相同技术手段的处理，对其技术路线和原理不再进行相同阐述，文中给出了对应信息的结果及其结果分析。

二、问题假设

1.忽略X射线的体积和形状，以质点处理。

2.忽略X射线的衍射，将其视为直线传播。

3.假设X射线由单一能量 (波长相同)的光子组成。

4.假设忽略探测器所探测到的一次或多次散射光子。

?

其中，λ1为通过小圆测得的探测器单元之间的距离，d为小圆的直径，n1为小圆所占探测器单元的个数。椭圆与小圆之间的距离为x=26mm，椭圆的长半轴为a=40mm，椭圆的短半轴为b=15mm。

三、模型的建立与求解

(一)问题一模型的建立与求解

CT系统的X射线是平行等距光线，整个发射－接受系统绕一固定中心逆时针方向旋转180次，在每一个X射线方向，具有512个等距单元探测器。

运用Excel将附件2中的数据进行处理，将大于0的数据单元格与等于0的数据单元格进行不同颜色区分标识，找出数据中的最大值、最小值、近似值。

得出结论1:基于透射原理分析可知，即探测器的接受信息数据最大值位置是X射线穿过椭圆的长轴处。

图像数据处理:将小圆数据和椭圆的数据非重叠区域，进行分离。

1.探测器单元之间距离问题求解

进行数据预处理，将单元格列宽调整为与行高相同的值，得到512×180个点，不难发现探测器逆时针共旋转了180°，起始方向由下斜向左上方向。由于正方形托盘上放置两个均匀固体介质，一个为椭圆，一个为半径为4mm的圆，由该模板的接收信息可以分析得到小圆的投影在旋转180次均为直径，用软件计算出共有29个接收点，共有28个探测器单元

其中，n2为椭圆边界与小圆边界之间最短距离情况下所占的探测器单元个数;n3为椭圆长轴所占的探测器单元个数;n4为椭圆短轴所占的探测器单元个数。求出ni的个数，计算出λi，然后求取平均值珔λ。珔λ即为探测器单元之间的距离。

计算得出各情况下所占的探测器单元个数:

n2=93;n3=288;n4=107

可以求解出各情况下的探测器单元之间的距离:

求取λi的平均值:

将各组数据代入上式得:

珔λ=0.2808mm

(二)问题二模型的建立与求解

1.模型建立

图1 滤波反投影的具体步骤

(1)滤波反投影算法:

Step1:对某 θ角度下的投影 pθ(t)作一维傅里叶变换，记为 Sθ(w);

Step2:对S0(w)乘以一维傅里叶权重因子ρ;

Step3:对第二步结果作一维傅里叶反变换，记作Qθ(t);

Step4:将0°～180°的所有修正过的投影函数Qθ(t)作为直接反投影计算，得到断层图像f(x，y)。

与滤波函数LRL(ρ)相对应的卷积函数lRL(R)为:

R－L滤波函数形式简单、实用，用它重建图像清晰。

2.模型求解

CT系统得到某未知介质的接收信息，由问题一的标定参数分析，即解决CT图像重建问题，问题分支为:

(1)物体在正方形托盘中的几何位置

(2)物体的几何形状

(3)物体的吸收率

几何位置、几何形状的求解。基于问题一的参数标定，建立滤波反投影重建模型，解决问题二中 (1)(2)问题分支，步骤如下:

Step1:加载数据，取投影

原始附件3数据给出了探测器在180个方向上的接收信息。因为在每个方向上，探测器接收的数据，即代表在这个方向上未知介质的投影信息。

Step2:滤波

使用Matlab里的Iradon函数中‘Ram－Lak’滤波器，对数据进行滤波运算处理。

Step3:坐标系变换

原始坐标系是以正方形托盘中心 (即问题一中的椭圆中心)为原点，现以CT系统旋转中心为原点，需将原始坐标系进行平移变换。

Step4:修正原始数据

填充探测器接收信息数据，为了更正确一条X射线在180个方向的吸收率与投影关系，将180个方向上的512组探测器单元填充至712组，填充数据用0来填充。即在上下分别增加一组100行180列的0矩阵。

Step5:Radon逆变化重建图像

将修正后的数据，在180个方向上进行Radon变化，同时改变坐标系的单位，其中

分别为X，Y轴的坐标系取值范围，假设在180个方向上，每次改变相同角度1度，进行反投影重建。并使用Matlab将正方形托盘模拟绘制在坐标系中。

Step6:针对几何位置和几何形状得出结果

从重建后的图像结合相关资料分析可知，未知介质的图像的几何形状是人类头部。即CT系统对人头部的扫描结果。其几何位置关系对应实际坐标系如图13所示。

(三)问题三模型求解

针对问题三的要求，提出了基于问题二的模型建立与求解的方法。

结果分析，提取未知介质的几何位置，几何形状，吸收率等相关信息进行定性定量分析:

A.几何位置:由图15可知反投影重建后的图像，分布在正方形托盘上的不均匀位置处，定性的分析，具体位置坐标结果结合附件problem3。

B.几何形状:图形分布较为涣散，是一种分布无规则的网织状图形

［1］郭立倩.CT系统标定与有限角度CT重建方法的研究 ［D］.大连理工大学，2016.

［2］吴孟达.数学建模教程 ［M］.高等教育出版社，2011.

［3］骆岩红.CT图像重建滤波反投影算法中指数滤波器的研究［J］.计算机科学，2014，41(S1):220－223.［2017－09－17］.

［4］张俊，闫镔，陆利忠，李磊，张峰.一种锥束CT系统几何参数标定方法［J］.核电子学与探测技术，2014，34(02):235－238+242.［2017－09－17］.

［5］黄亚，张祥志，祝江威，刘海岗，邰仁忠.CT断层图像重建的新滤波函数 ［J］.核电子学与探测技术，2012，32(12):1388－1393.［2017－09－17］.

荣亚亚 (1996.5－)，女，汉族，河北省石家庄市，华北理工大学经济学院经济统计学专业;郭子健 (1998.3－)，男，汉族，湖南省冷水江市，华北理工大学理学院，应用统计学专业;吉督(1997.10－)男，汉族，湖南省长沙市浏阳，华北理工大学机械工程学院，机械设计制造及其自动化;赵旭东 (1996.11－)，男，汉族，河北省石家庄市。华北理工大学建筑工程学院，物流工程专业。