论文部分内容阅读
[摘 要]本文針对已知CT系统基本信息，通过数学分析、matlab和Excel运算、FBP反演公式和iradon函数对其几何参数进行标定，并进行精度、稳定性分析及改进探讨 。CT系统可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息，但是CT系统安装时往往存在误差，从而影响成像质量。
[关键词]CT系统、FBP反演、iradon函数、增益
中图分类号：G717 文献标识码：A 文章编号：1009-914X（2018）48-0023-02
一、模型假设
假设正方形托盘绝对水平，托盘在转动过程中两个均匀固体始终保持绝对静止；
假设X射线发射器和探测器的旋转为等角度旋转；
假设旋转中心在探测器组所在线段的垂直平分线上；
假设每个探测器为一个点，不考虑探测器本身宽度；
假设探测器排列顺序为沿顺时针方向从1-512顺序排列；
假设两个均匀固体介质理想化，忽视其厚度，将其分别视为二维平面中的标准椭圆和圆；
假设模板接收信息的第一列数据为CT系统经第一次旋转前测量并处理的信息；
假设问题二反演出的几何图形中心和正方形托盘中心重合。
二、符号说明
d 探测器间距
l 圆形介质的直径
m接收窄带信号的探测器个数
p最宽色带中点向y轴的投影长度
q最窄色带中点向x轴的投影长度 αCT系统每次转过的角度（步进角度）
βx射线初始出射角
θ切线a逆时针方向与x正半轴夹角
切线b逆时针方向与x正半轴夹角
注：未说明符号在文中用到时注明
三、构建模型
3.1图形分析
首先在模板示意图以正方形中心（即椭圆介质的圆心）为原点、水平中轴线为x轴建立直角坐标系，如图1所示；对附件二做底色填充处理，将数据按升序（为0的单元格除外）由红黄双色阶过渡填充，得出探测器所接收的射线能量分布图2。分析图2可知：
1）由图对应模板形状可知，较宽的着色部分为x射线穿过椭圆介质后探测器的接收信息；似窄带形着色部分为x射线穿过圆介质后探测器的接收信息，因圆为中心对称图形，故窄带宽度近似相等。
2）有色部分最宽处是X射线平行x轴入射，最窄处为沿y轴平行入射。
3）探测器在初始位置和180次扫描后两次接收的射线能量分布关于整体图形中心近似对称，结合模板本体为中心对称图形，可猜想扫描180次近似旋转180度。此猜想会在接下来作详细证明。
4）a点为x射线沿图1 a路线射入时的临界点，从此点往后两着色带将逐渐重叠直至完全覆盖；b点为x射线沿图1.b路线射入时的临界点，从此点往后两着色带将开始分离直至完全分离。
3.2探测器间距求解
在已确定窄带着色区为x射线穿过圆形介质后的接收信息前提下，探测器间距d可由圆形介质直径l与接收窄带信号的探测器个数m作商求得，即
（1）
分析图2可得，接收窄带信号的探测器个数m为29，圆形介质直径l由模板几何信息可得为8mm，故求得探测器间距
3.3旋转中心求解
旋转中心是通过对图2较宽的着色区最窄和最宽两条色带所对应的探测器位置分析得到的（即x射线水平和垂直射向椭圆介质时）。对于最窄和最宽两条色带，通过Excel函数分别找到最宽和最窄色带的中点p，q，并记录相对应的行数，即第几个探测器。结合假设2旋转中心始终在探测器第256个单元和第257个单元之间的平分点且过该中心平行于探测器方向的直线上，和探测器间距，即可求得旋转中心到椭圆介质长轴和短轴的投影长度，从而求得旋转中心在已构建直角坐标系中的坐标值。
Excel函数为：512-（COUNTIF（R[-512]C：R[-1]C，"0"）），运算结果如图3：
对图3第513行数据筛选标记，筛出最小值标记为蓝色（如图4），最大值标记为红色（如图5），所在色带的中点所对应的探测器单元个数（若标记数据个数大于等于2则取中间列，如图5阴影所示）。
经分析计算得，最宽色带中点对应第235个探测器单元，最窄色带对应第224个探测器单元。结合探测器间距计算可得最宽色带中点向y轴的投影长度p和最窄色带中点向x轴的投影长度q如下：
结合假设4探测器排列规则，可确定旋转中心在第二象限，坐标为（-8.69085，+5.65595）。具体位置如图6所示，红星标注点即为旋转中心。
3.4CT系统x射线转过的180个方向求解
结合假设1和假设6，欲求x射线转过的180个方向，需先求出CT系统x射线第一次的出射位置和每次转过的角度。
3.4.1 CT系统每次转过的角度（步进角度）α
由3.2.1的第四点分析已知图1的两条切线分别对应图2 a、b两点。构建matlab算法，算法思想为首先结合模板几何信息得出椭圆介质的切线方程，利用圆与直线相切的条件为点到直线距离等于半径的性质，将圆介质半径 代入方程，即可得与椭圆介质和圆介质均相切的a、b两条直线，利用两条直线斜率即可得出a、b直线逆时针方向与x正半轴所夹角度 和 。对应附件二确定a、b对应次数分别为第14次和109次，差值为95次，输入算法中，运行即可求得步进角度 为 。
关键代码如下：
syms x y；
[x，y]=solve（'（x/5-1）^2=16*（（x/225）^2+（y/1600）^2）'，'x^2/225+y^2/1600=1'，'x'，'y'）；x=max（x）； y=sqrt（（1-x^2/225）*1600）；
beta=2*rad2deg（double（atan（1600*x/（225*y））））；
运行结果如下：
angle =
1.0005
由4.2.1中分析第一点可知，最宽着色带所对应的是x射线平行于是x轴入射，有两种可能性，故在假设4对探测器设定前提下，入射方向为与x轴正方向呈 。因系统逆时针旋转，结合假设6，只需确定最宽着色带所对应的扫描次数x，即可得到x射线第一次出射位置与x轴正半轴夹角即初始出射角β，公式如下： （2）
查附件二可知，扫描次数x=61，则可得初始出射角β= ，可以确定x射线初次出射方向为从第四象限 第二象限，与x轴正方向所夹角度为 。
四、未知物质反演推导
4.1 几何介质图形的反演求解
通过问题一已知此CT系统的旋转中心、步进角度和探测器间距等基本信息。在此基础上，欲通过用上述CT系统得到的某未知介质的接收信息，反演确定该未知介质在正方形托盘上的相关信息。常见反演公式有平面束投影变换及其反演、FBP反演公式、BPF反演公式等，我们选择FBP反演公式来计算。
FBP反演公式即为滤波反投影重建公式，是在实际CT图形重建中最为常用的一种公式 （3）
利用此公式和iradon函数，实现重建图像的过程如下：
把投影矩阵R转换到频域，生成fft（R）；
fft（R）和滤波函数H相乘，得到滤波后的频域投影矩阵fft（R）*H；
把fft（R）*H转换到空域，得到空域中的滤波后的投影矩阵 = ifft（fft（R）*H）；
插值后，得到处理好的投影矩阵 ；
直接反投影得到重建图像矩阵。
关键代码如下I=iradon（R，theta，interp，filter，frequency_scaling，output_size）
具体参数介绍：
R是投影矩阵；
theta是一个包含所有扫描角度的向量，这时每两个相邻角度间隔角度间隔相等；也可以是一个标量值，等于相邻两个扫描角度的间隔；
interp插值函数；
filter是滤波函数；
frequency_scaling是一个标量值，取值范围[0，1]，通过缩减滤波函数的频率修改滤波函数。Output_size是一个标量，用来规定重建图像的行数和列数。
将附件三数据代入matlab中运行，得图7：
观察图7可知所求几何图形占图像的部分，为了得到256*256的像素矩阵，故通过改进matlab算法，截取以介质几何图形为中心的362*362的像素矩阵，并将其缩放、旋转成256*256的像素矩阵，运行结果如图8所示：
详细代码请见附录一。
五、模型的检验与评价
5.1验证
5.1.1旋转角度的验证
在excel中通过对探测器180次扫描数据的分析，可以得出：第151次扫描与第61次扫描相差90度，而这两次之间相差90次扫描，故得出平均一次扫描旋转角度为1度，与求得角度1.0005度几乎吻合，验证成功。
5.1.2 FBP模型的验证
用第一问求得的标定参数通过iradon函数重建图形，得到如右下图形：
将该图与附件1所给数据拟合得到的图形（左上）比较，可以看出图形几乎重合，因此iradon函数重建模板的算法验证成功。
5.2 评价
在CT成像系统的几何参数获取过程中，误差的存在不可避免。误差的主要来源分为两种：机械误差精度和标定模板投影图像的像素坐标获取精度，后者远大于与前者。
对于旋转中心等标定参数来说，数学分析相较于正弦图拟合在精度方面有一定的误差。
iradon函数滤波反投影重建模型相较于代数重建算法等其他离散型CT模型在滤波方面有更好的拟合性。
求解模板相应点的吸收率时，可以采用拟合的方法，利用周围的值拟合出该点的值，误差会相对减小。
将362*362像素矩阵缩放为256*256像素矩阵时，除了使用线性插值法外，还可以使用最近邻域法缩放，这样可以进一步改善图像的失真度，最大程度还原图像。
在使用iradon函數重建图形时，增大重建投影矩阵的容量可以防止矩阵顶点处数据还原失真，有效提高精度。
参考文献
[1]司守奎，孙兆亮.数学建模算法与应用（2版）.[M].北京：国防工业出版社，2016.
[2]（美）米尔斯切特著；刘来福，黄海洋，杨淳译.数学建模方法与分析.[M].北京：机械工业出版社，2014.
[3]姜启源，谢金星，叶俊.数学模型（4版）.[M].北京：高等教育出版社，2011.
[4]胡良剑.MATLAB数学实验.[M].北京：高等教育出版社，2011.
[5]闫镔，李磊.CT图像重建算法.[M].北京：科学出版社，2014.
[6]张朋，张慧滔，赵云松.X射线CT成像的数学模型及其有关问题[D].北京：首都师范大学数学科学学院 检测成像工程研究中心.2012：5-6.
[关键词]CT系统、FBP反演、iradon函数、增益
中图分类号：G717 文献标识码：A 文章编号：1009-914X（2018）48-0023-02
一、模型假设
假设正方形托盘绝对水平，托盘在转动过程中两个均匀固体始终保持绝对静止；
假设X射线发射器和探测器的旋转为等角度旋转；
假设旋转中心在探测器组所在线段的垂直平分线上；
假设每个探测器为一个点，不考虑探测器本身宽度；
假设探测器排列顺序为沿顺时针方向从1-512顺序排列；
假设两个均匀固体介质理想化，忽视其厚度，将其分别视为二维平面中的标准椭圆和圆；
假设模板接收信息的第一列数据为CT系统经第一次旋转前测量并处理的信息；
假设问题二反演出的几何图形中心和正方形托盘中心重合。
二、符号说明
d 探测器间距
l 圆形介质的直径
m接收窄带信号的探测器个数
p最宽色带中点向y轴的投影长度
q最窄色带中点向x轴的投影长度 αCT系统每次转过的角度（步进角度）
βx射线初始出射角
θ切线a逆时针方向与x正半轴夹角
切线b逆时针方向与x正半轴夹角
注：未说明符号在文中用到时注明
三、构建模型
3.1图形分析
首先在模板示意图以正方形中心（即椭圆介质的圆心）为原点、水平中轴线为x轴建立直角坐标系，如图1所示；对附件二做底色填充处理，将数据按升序（为0的单元格除外）由红黄双色阶过渡填充，得出探测器所接收的射线能量分布图2。分析图2可知：
1）由图对应模板形状可知，较宽的着色部分为x射线穿过椭圆介质后探测器的接收信息；似窄带形着色部分为x射线穿过圆介质后探测器的接收信息，因圆为中心对称图形，故窄带宽度近似相等。
2）有色部分最宽处是X射线平行x轴入射，最窄处为沿y轴平行入射。
3）探测器在初始位置和180次扫描后两次接收的射线能量分布关于整体图形中心近似对称，结合模板本体为中心对称图形，可猜想扫描180次近似旋转180度。此猜想会在接下来作详细证明。
4）a点为x射线沿图1 a路线射入时的临界点，从此点往后两着色带将逐渐重叠直至完全覆盖；b点为x射线沿图1.b路线射入时的临界点，从此点往后两着色带将开始分离直至完全分离。
3.2探测器间距求解
在已确定窄带着色区为x射线穿过圆形介质后的接收信息前提下，探测器间距d可由圆形介质直径l与接收窄带信号的探测器个数m作商求得，即
（1）
分析图2可得，接收窄带信号的探测器个数m为29，圆形介质直径l由模板几何信息可得为8mm，故求得探测器间距
3.3旋转中心求解
旋转中心是通过对图2较宽的着色区最窄和最宽两条色带所对应的探测器位置分析得到的（即x射线水平和垂直射向椭圆介质时）。对于最窄和最宽两条色带，通过Excel函数分别找到最宽和最窄色带的中点p，q，并记录相对应的行数，即第几个探测器。结合假设2旋转中心始终在探测器第256个单元和第257个单元之间的平分点且过该中心平行于探测器方向的直线上，和探测器间距，即可求得旋转中心到椭圆介质长轴和短轴的投影长度，从而求得旋转中心在已构建直角坐标系中的坐标值。
Excel函数为：512-（COUNTIF（R[-512]C：R[-1]C，"0"）），运算结果如图3：
对图3第513行数据筛选标记，筛出最小值标记为蓝色（如图4），最大值标记为红色（如图5），所在色带的中点所对应的探测器单元个数（若标记数据个数大于等于2则取中间列，如图5阴影所示）。
经分析计算得，最宽色带中点对应第235个探测器单元，最窄色带对应第224个探测器单元。结合探测器间距计算可得最宽色带中点向y轴的投影长度p和最窄色带中点向x轴的投影长度q如下：
结合假设4探测器排列规则，可确定旋转中心在第二象限，坐标为（-8.69085，+5.65595）。具体位置如图6所示，红星标注点即为旋转中心。
3.4CT系统x射线转过的180个方向求解
结合假设1和假设6，欲求x射线转过的180个方向，需先求出CT系统x射线第一次的出射位置和每次转过的角度。
3.4.1 CT系统每次转过的角度（步进角度）α
由3.2.1的第四点分析已知图1的两条切线分别对应图2 a、b两点。构建matlab算法，算法思想为首先结合模板几何信息得出椭圆介质的切线方程，利用圆与直线相切的条件为点到直线距离等于半径的性质，将圆介质半径 代入方程，即可得与椭圆介质和圆介质均相切的a、b两条直线，利用两条直线斜率即可得出a、b直线逆时针方向与x正半轴所夹角度 和 。对应附件二确定a、b对应次数分别为第14次和109次，差值为95次，输入算法中，运行即可求得步进角度 为 。
关键代码如下：
syms x y；
[x，y]=solve（'（x/5-1）^2=16*（（x/225）^2+（y/1600）^2）'，'x^2/225+y^2/1600=1'，'x'，'y'）；x=max（x）； y=sqrt（（1-x^2/225）*1600）；
beta=2*rad2deg（double（atan（1600*x/（225*y））））；
运行结果如下：
angle =
1.0005
由4.2.1中分析第一点可知，最宽着色带所对应的是x射线平行于是x轴入射，有两种可能性，故在假设4对探测器设定前提下，入射方向为与x轴正方向呈 。因系统逆时针旋转，结合假设6，只需确定最宽着色带所对应的扫描次数x，即可得到x射线第一次出射位置与x轴正半轴夹角即初始出射角β，公式如下： （2）
查附件二可知，扫描次数x=61，则可得初始出射角β= ，可以确定x射线初次出射方向为从第四象限 第二象限，与x轴正方向所夹角度为 。
四、未知物质反演推导
4.1 几何介质图形的反演求解
通过问题一已知此CT系统的旋转中心、步进角度和探测器间距等基本信息。在此基础上，欲通过用上述CT系统得到的某未知介质的接收信息，反演确定该未知介质在正方形托盘上的相关信息。常见反演公式有平面束投影变换及其反演、FBP反演公式、BPF反演公式等，我们选择FBP反演公式来计算。
FBP反演公式即为滤波反投影重建公式，是在实际CT图形重建中最为常用的一种公式 （3）
利用此公式和iradon函数，实现重建图像的过程如下：
把投影矩阵R转换到频域，生成fft（R）；
fft（R）和滤波函数H相乘，得到滤波后的频域投影矩阵fft（R）*H；
把fft（R）*H转换到空域，得到空域中的滤波后的投影矩阵 = ifft（fft（R）*H）；
插值后，得到处理好的投影矩阵 ；
直接反投影得到重建图像矩阵。
关键代码如下I=iradon（R，theta，interp，filter，frequency_scaling，output_size）
具体参数介绍：
R是投影矩阵；
theta是一个包含所有扫描角度的向量，这时每两个相邻角度间隔角度间隔相等；也可以是一个标量值，等于相邻两个扫描角度的间隔；
interp插值函数；
filter是滤波函数；
frequency_scaling是一个标量值，取值范围[0，1]，通过缩减滤波函数的频率修改滤波函数。Output_size是一个标量，用来规定重建图像的行数和列数。
将附件三数据代入matlab中运行，得图7：
观察图7可知所求几何图形占图像的部分，为了得到256*256的像素矩阵，故通过改进matlab算法，截取以介质几何图形为中心的362*362的像素矩阵，并将其缩放、旋转成256*256的像素矩阵，运行结果如图8所示：
详细代码请见附录一。
五、模型的检验与评价
5.1验证
5.1.1旋转角度的验证
在excel中通过对探测器180次扫描数据的分析，可以得出：第151次扫描与第61次扫描相差90度，而这两次之间相差90次扫描，故得出平均一次扫描旋转角度为1度，与求得角度1.0005度几乎吻合，验证成功。
5.1.2 FBP模型的验证
用第一问求得的标定参数通过iradon函数重建图形，得到如右下图形：
将该图与附件1所给数据拟合得到的图形（左上）比较，可以看出图形几乎重合，因此iradon函数重建模板的算法验证成功。
5.2 评价
在CT成像系统的几何参数获取过程中，误差的存在不可避免。误差的主要来源分为两种：机械误差精度和标定模板投影图像的像素坐标获取精度，后者远大于与前者。
对于旋转中心等标定参数来说，数学分析相较于正弦图拟合在精度方面有一定的误差。
iradon函数滤波反投影重建模型相较于代数重建算法等其他离散型CT模型在滤波方面有更好的拟合性。
求解模板相应点的吸收率时，可以采用拟合的方法，利用周围的值拟合出该点的值，误差会相对减小。
将362*362像素矩阵缩放为256*256像素矩阵时，除了使用线性插值法外，还可以使用最近邻域法缩放，这样可以进一步改善图像的失真度，最大程度还原图像。
在使用iradon函數重建图形时，增大重建投影矩阵的容量可以防止矩阵顶点处数据还原失真，有效提高精度。
参考文献
[1]司守奎，孙兆亮.数学建模算法与应用（2版）.[M].北京：国防工业出版社，2016.
[2]（美）米尔斯切特著；刘来福，黄海洋，杨淳译.数学建模方法与分析.[M].北京：机械工业出版社，2014.
[3]姜启源，谢金星，叶俊.数学模型（4版）.[M].北京：高等教育出版社，2011.
[4]胡良剑.MATLAB数学实验.[M].北京：高等教育出版社，2011.
[5]闫镔，李磊.CT图像重建算法.[M].北京：科学出版社，2014.
[6]张朋，张慧滔，赵云松.X射线CT成像的数学模型及其有关问题[D].北京：首都师范大学数学科学学院 检测成像工程研究中心.2012：5-6.