数学建模大赛CT系统参数的标定与成像问题.docx
《数学建模大赛CT系统参数的标定与成像问题.docx》由会员分享，可在线阅读，更多相关《数学建模大赛CT系统参数的标定与成像问题.docx（49页珍藏版）》请在冰点文库上搜索。
数学建模大赛CT系统参数的标定与成像问题
CT系统参数标定及成像
芋
本文主要研究CT系统参数的标定与成像问题，通过建立模型利用结构样品标定CT系统的参数，并对未知结构样品的图像进行重建。
针对问题一，先对圆分析，无论'射线的方向如何，所得到的投影长度总是等于圆的直径，由此关系可求得相邻探测器单元的距离为。
再利用几何知识求出圆和椭圆在某个投影方向上的投影长度与该方向角度的关系式，并将其有较大跳跃或异常的点去除，利用傅里叶变换拟合出一个函数，再将附件2中180个角度对应的数据与拟合函数在同一角度得出的数据做差，用最小二乘法拟合出一个更加精确的函数，再将附件2中的数据代入进去求出180个角度，其中前三个以y轴正方向的夹角分别是29.0661°,30.5702°,31.5755°。
以托盘左下角为原点，水平方向为横坐标，竖直方向为纵坐标建立直角坐标系，同样利用最小二乘法得出旋转中心为（40.7543,56.2520）。
针对问题二，首先找到投影坐标与原图像坐标之间的转化关系，再利用滤波反投影法得出图像的衰减系数分布函数，即所有经过滤波后的射线投影在0-180度范围的衰减系数的叠加。
考虑到CT系统是逆时针旋转180次，每一次对应一个角度，x射线对物体进行平行束扫描，得到是180个离散投影数据，于是将衰减系数分布函数离散化，将180个角度得到的投影数据叠加，得到图像的像素点，从而重建出图像，并算出了附件四给出的10个位置的吸收率分别为:
、、、、、、、、、0.0307o
针对问题三，利用问题二中的模型重建出图像并解出10个位置的吸收率分别为：
、、、、、、、、、0.0774o
针对问题四，关于参数标定的精度与稳定性问题，利用CT系统的结构以及成像原理可知，x射线对物体进行平行束扫描获得了180个离散数据，由于CT系统的探测器之间具有一定的间隔，所以物体的边界可能有一局部无法被收集到信息，从而造成投影长度缩小的情况，导致所求探测器的单元间隔的增大，并且造成旋转中心与x射线旋转角的误差。
通过分析椭圆投影到探测器上长度所对应的误差分析，得到投影长度越长，模板以下图形面积越小，所对应的误差就越小，即可得出一个精度较高的新模板。
关键词：
最小二乘法；滤波反射法；衰变系数；傅里叶变换;
一、问题的重述
问题的背景：
计算机层析成像（ComputedTomography,CT）是通过对物体进行不同角度的射线投影测量而获取物体横截面信息的成像技术，涉及多个学科领域。
CT技术不但给诊断医学带来革命性的影响.还成功地应用于无损检测、产品反求和材料组织分析等工业领域，其实质是由扫描所得到的投影数据反求出成像平面上每个点的衰减系数值。
有一种典型的二维CT系统，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。
X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。
对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信，息、O
问题的目的：
CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于结构的样品（称为模板）标定CT系统的参数,并据此对未知结构的样品进行成像。
所要解决的问题：
问题一：
在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息已经给定，相应的数据文件见附件1,其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。
对应于该模板的接收信息见附件2。
请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。
问题二：
附件3是利用上述CT系统得到的某未知介质的接收信息。
利用⑴中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。
另外，请具体给出所给的10个位置处的吸收率，相应的数据文件见附件4。
问题三：
附件5是利用上述CT系统得到的另一个未知介质的接收信息。
利用
（1）中得到的标定参数，给出该未知介质的相关信息。
另外，请具体给出所给出的10个位置处的吸收率。
问题四：
分析
（1）中参数标定的精度和稳定性。
在此根底上自行设计新模板、建立对应的标定模型，以改良标定精度和稳定性，并说明理由。
二、问题的分析
本文主要研究CT系统两方面的问题，一是利用结构的样品和它的投影
信息得到CT系统的参数，二是用己知的标定CT系统的参数来对未知结构的样品进行成像。
2.1问题一的分析：
该问题要求利用附件1和附件2的信息来确定出CT系统的旋转中心、相邻探测单元之间的距离以及X射线每次旋转的方向。
针对这个问题，考虑到圆的图形比拟规那么且投影也是规那么的，首先对圆单独进行考虑，对于任何一个方向的X射线，圆对应的投影长度都等于圆的直径长，分析可知圆的投影长度与接收到该圆信息的探测器单元的个数有关，由于每个探测器单元是等间距，由此计算出相邻探测器单元间的距离。
附件2中180列数据即是对应的每个x射线方向的探测器上所接受到的信息。
因为附件2中的每一列数据表示的都表示每次转角所对应的投影长度，用几何方法求解出的函数的那些突变值会对函数造成较大的影响。
理论上用附件2的数据拟合和几何计算出的射线转角与投影长度的函数关系应该是吻合的，但由于去除特殊点所拟合的函数也具有一定的误差，用数据拟合的函数会与其有一定的差值，为了使数据所拟合的函数与几何方法拟合的函数更加接近，可以最小二乘法进行拟合。
2.2问题二的分析：
该问题实际上是CT图像重建问题，即利用CT系统扫描得到的投影数据解出投影像平面上各像素点的衰减系数，然后重建CT图像。
由中心切片定理知道密度函数在某一方向上的投影函数的一维傅立叶变换⑵函数是原密度函数的二维傅立叶变换函数平面上沿同一方向且过原点的直线上的值。
首先找到投影坐标与原图形坐标之间的转化关系，然后求出某一角度下，在投影坐标系下X射线方向上衰减系数的线积分，易知投影方向不同或x射线离旋转中心的距离不同，x射线穿过模板的厚度那么不同，那么衰减系数的线积分也会不同，所以衰减系数的线积分是与x射线的角度。
以及投影射线到旋转中心的距离p有关的一个函数。
然后对其进行一维傅里叶变换，并利用坐标变换，转化为在原图像坐标系下的衰减系数的线积分函数，再对其进行二维的傅里叶变换⑵，并进行二维傅里叶逆变换，然后经过滤波函数得到滤波后的投影像素点的衰减系数的一个分布函数，该分布函数对应的是经过滤波后的射线投影在。
〜180度范围的衰减系数的叠加。
但实际在CT扫描时，得到是180个离散投影数据，而不是连续的角度变化，于是将衰减系数分布函数离散化，将180个角度得到的投影数据叠加，得到图像的像素点，从而重建出图像。
2.3问题三的分析：
问题三的求解过程与问题二的相同，都是利用滤波反射模型进行求解，从重建出图像和10个位置的吸收率。
2.4问题四的分析：
要分析参数的精度与稳定性，首先找到造成参数误差的原因，即分别对问题一中所求的CT系统的旋转中心的位置、探测器单元间隔以及x射线旋转方向进行误差分析。
要进行误差分析，可以从CT系统的结构以及成像原理入手，由于探测器之间具有一定间隔以及投影得到的是离散数据，分析可得，当x射线进行平行束扫描时，被投影物体的两侧的数据可能有一局部无法被接收器接受，从而导致投影的长度并不会等于对应原图像的长度，而是会比原像小一点，从而造成参数误差，再从这个误差来源入手，分析问题一的参数的精度问题。
三、模型的假设
1.假设X射线不会因外界其它因数的干扰而衰落；
2.假设每个探测器单元的长度都是相等的；
3.假设每束］射线都是一样的；
4.x射线的初始位置的方向是朝左上的；
四、符号的说明
符号
符号说明
圆的直径的大小
第，次旋转后，'射线射向椭圆在探测器上对应非“0”值
的探测器单元的个数
第，次旋转后，'射线射向圆在探测器上对应非“0”值的
探测器单元的个数
AZ
每个探测器的长度
d总
椭圆和圆的投影总长度
五、模型的建立与求解
5.1对于问题一的求解
5.1.1问题一的解题步骤：
（1）根据附件2模板接受的球的信息来确定相邻两探测器单元之间的距离；
（2）随x射线旋转角度的变化，利用几何知识列出探测器上由椭圆和圆吸收x射线所导致的非“0”值探测器单元的个数与x射线方向的函数关系式。
再用附件2中每一列非“0”值的探测器单元的个数，用最小二乘法来拟合出x射线每个角度与对应的非“0”值探测器单元的个数的函数再与几何求解出的函数进行比拟，从而确定x射线方向与探测器上由椭圆和圆吸收'射线所导致的非“0”值探测器单元的个数的函数关系式；
（3）利用最小二乘法求解旋转中心的位置。
模型的建立与求解
5.1.2.1x射线入射方向模型
（一）相邻两探测器单元之间距离确实定
附件1中的数据中只有“0”和“1”的两种值，将其导入到MATLAB中可以发现，“1”值对应的为椭圆和圆的投影数据，“0”值那么代表正方形托盘的其他的地方的投影数据，由此可以判断出问题给出的“吸收率”是指固体介质的对x射线的吸收率，并且“吸收率”为1（即可以认为投射到固体介质上的所有x射线都被吸收了）。
由于］射线是平行照射到圆上，因此无论射线从任何方向照射圆，对应角度的探测器上的非“0”值的探测器的个数都是一个定值并且这些探测器单元的总长度都应该为球的直径的长度。
附件2中的每列数据是对应每次旋转后的非“0”探测器单元的个数，理论上对于每列圆所对应于非“0”探测器单元的个数应该相等。
但是附件2中每一列圆对应于非“0”探测器的个数不相等，导致这种现象的原因可能是有些时候探测器感应出现问题、x射线在传播过程中发生了折射等因素导致的，如果用非“0”探测器的个数用最大或者最小值来代替会出现取的最大或最小值是由意外情况导致的，因此用平均值来进行计算得到每个探测器单元的长度为：
△/=斗=0.2775
A
i=l
其中N为圆在探测器上投影的列数。
（二）x射线的180个方向确实定
1.几何方法确定180个方向与投影总长度的关系
根据图一所示的附件2中x射线在探测器上的投影图（竖直方向的白色长度集合表示X射线以不同入射角照射椭圆或圆投射在探测器上接受到非“0”值的所有探测器单元的总长度集合），为了更好的表述，定义椭圆投影长度为椭圆在某个投影方向上所对应的图一中的一段竖直的白色长度，圆投影长度为圆在某个投影方向上所对应的图一中所对应的一段竖直的白色长度。
图一附件2中x射线在探测器上的投影图
以椭圆中心为原点，水平向圆心的方向为x轴的正方向，竖直向上的方向为y轴的正方向，建立直角坐标系，初始位置的］射线入射方向向左上，射线的方向与y轴正方向的夹角为。
角。
1由初始位置（图像最左边）到A点椭圆的投影长度
根据图一可以发现由图像最左边到A点的这段过程中在同一个竖直方向上有两段白色长度（即椭圆与圆的投影长度没有重叠），白色的长度即可认为是椭圆和圆投影长度的和。
0
L:
Y=kX+b0
x轴
图二X射线只与椭圆相切图
当入射的X射线只与椭圆相切（如图二所示）时假设入射的X射线乙的方程为:
y=kx+b0；
当入射的工射线与椭圆相切时有：
y=kx+b0
F——=1
U52402
联立上述的两个方程组有：
（1600+225^2）%2+450^x+225Z?
02—360000=0；
因为入射的X射线与椭圆相切，所以有：
△=（450002—4（1600+225尸）（225/?
。
2-360000）=0；
得到:
225炉―膏+1600=0；
b：
=225好+1600；
从而求解出入射的尤射线与椭圆相切的两个切点分别为:
/.-225k1600、
（玉，y1）=（/y，I；
丁1600+225好J1600+225#2
、/.225-1600、
（心，yQ=（I==,I=）；
~一J1600+225好J1600+225仃
那么这两条切线所对应的方程分别为：
乂=kx+J225比之+1600;
%=奴—.225妃+1600；
那么这两条切线（椭圆投影长度）的距离为：
-7225k2+1600-7225k2+1600
d=（其中R=tan（90。
+。
））；
W+1
代入得：
2^225tan2（90°+6>）+1600
d息=I.+d（）；
Vtan2（90°+6>）+1
根据两条切线的距离再结合附件2中给出的第一列数据中投影的长度，来求解出射线的初始入射角度大约为29.25。
。
2由A点到C点位置椭圆的投影长度
根据图一可以发现竖直方向的两端白色带合为了一条，说明了椭圆和圆的投影长度有重叠，即图二所示，圆只有一局部吸收了x射线，另一局部的］射线被椭圆吸收了所以导致椭圆和圆的投影长度有重合。
y轴
图三射线与圆相交于椭圆相切
在图三椭圆的上半局部，设有直线A和椭圆相切，直线七2和圆相切，且两直线互相平行，两直线的斜率分别电，两直线与y轴正方向的夹角为劣。
直线％和直线人的方程分别为：
yFx+九；V=k\X+b?
;
根据直线A和椭圆相切得出以下方程：
y=
U52402
联立以上两个式子得：
△=225妒+1600-疽=0
根据直线&2和圆相切得出以下方程：
Jy=kxx+b2\x-45）2+v2=16
联立以上两个式子得：
△=—20XXZ?
—90—仇2+16=0
根据（3）和（4）式子可得：
b、=J225妃+1600；
-90#|+J8100#，—4（20XX上「16）
。
2=2；
解得：
k、=—1.0922或—
所以直线A和直线七2的距离（即圆的投影长度）为：
如-b]
2依+1
阻-b\2』225妒+1600+90k、—^SlOO^2-4（20XX^2-16）d=z
也2+1
2」225tan2（90°+^）+1600+90tan（90°+%）」8100tan2（90°+%）—4（20XX+tan（90°+^）-16）
2^tan2（90°+^）+l
椭圆的投影长度和①中的求法一致为:
（―1.0922-0.8248）
所以投影的总长度为：
dz=d+d[°
在椭圆的下半局部求解一直线和椭圆相切，另一直线和圆相切，且两直线互相平行，这两直线的斜率和两直线与y轴正方向的夹角的求解和椭圆的上半局部的求解方法一样，解出的直线斜率为心或，总的投影长度也是圆和椭圆的投影长度的和。
3C点到边界的最右端
由于C点到边界的最右端根据图一可以发现由图像最左边到A点的这段过程中在同一个竖直方向上有两段白色长度（即椭圆与圆的投影长度没有重叠），白色的长度即可认为是椭圆和圆投影长度的和为d总-d+dG在
（0.81481.0921）上成立。
综上所述：
2V225F+1600,r^=+4，
如+1
2^225好+1600
TF+i，
2】225站+1600+90佑一^8100/:
/-4（20XX/:
；-16）2a/225^2+1600
'二，
（其中A=tan（9O°+0））,d°=8mm,d总表示180个射线方向角与圆和椭圆的
投影总长度的关系）
对上述函数先将有异常的点去除，再用傅里叶五级变换进行对其函数进行拟合，得到以下图像：
图四傅里叶五级变换后的函数
2.用附件2数据确定180个方向与投影总长度的关系
对于附件2中所给出的512*180的数据，可以知道每一列即为对应的一个X射线转动角度，每一列数据中非“0”值的个数即为圆和椭圆的投影总长度（X射线被圆和椭圆吸收后在探测器上射线信号减落的探测器单元的总个数）。
为了描述X射线的入射角度的变化对圆和椭圆的投影总长度（一列上非“0”数值的个数）的变化关系，可以用每一列上的非“0”值的个数做为纵坐标，以列数（x射线的每次转动的方向）做为横坐标，用最小二乘法来拟合出一个最正确的函数来。
设寻表示第i列数，y,•表示的是第i列数中非“0”值的个数，每一列（X射线的每次转动的方向）所对应的（X,,北）那么构成了含有180个点的点集，并且二者满足关系叫=f'CJCiy=a1xm+a2xm~1+……+洋"0冲1。
现要根据附件2的数据来拟合出一个次数不超过〃2次的多项式g（x）二,根据上述条件
/=0
即求解：
/=o
所解出的函数即为每次转动X射线所对应的角度与投影总长度的函数关系式。
3.x射线方向与投影总长度关系的模型
上文分别应用了几何方法求解函数关系和最小二乘法拟合函数的方法来确定的］射线的方向与投影总长度之间的关系。
这两种方法求解出的函数应该是同一个函数因为都是表示X射线的方向与投影总长度之间的函数关系，只是这两种的方法在求解过程中会有一定的误差，可能导致两函数可能在相同的］射线的方向上的投影总长度之间的函数值会有所偏差。
因为用几何计算出的180个方向与投影长度的关系有一些异常的点，因此可以将那些异常的点去除对那些正常的点进行拟合，能拟合出函数为』总=乙（9）；数据是用最小二乘法拟合，拟合出函数为d总=g（。
）。
由于几何计算出的点有些是异常点，需要进行拟合，在拟合过程中会有很多不同的拟合方法，用数据拟合也是如此。
因为上述两种方法求解的函数是同一个函数，理论上值应该相等，为了找出一个能更加准确地表示x射线方向与投影总长度关系，在所有相同偏转角的情况下两种方法拟合的函数所对应的值的差值的平方应该尽肯能的小。
因此可以建立x射线方向与投影总长度关系的模型，模型如下所示：
180
而（Z（乙⑶-幺⑶）
i=l
用MATALAB对进行模型求解后，得到了一个较优的射线方向与投影总长度的函数关系：
图五优化后射线方向与投影总长度的函数
再将附件2中每一列非“0”值的个数分别代入上述方程从而解得这180个角度方向如表一（详细的数据见附录）示所示：
表一X射线的180次的方向
编号
1
2
3
178
179
180
与y轴正方向的夹角
29.0662
30.5702
31.5728
207.0412
208.0438
209.0465
旋转中心模型
/\
0
坐标票点
图六旋转图
较转中心
（W3）
如图六中所示，以Y。
为坐标原点，工射线射向的方向为水平方向的正方向
建立直角坐标系。
由旋转几何关系可以推断出直线FoK的长度可以用随角度变化的函数尸3（。
）来表示其长度：
气*=匕'=如0）;
设Y是一个常量点（行匕）与点（X39y?
两点连线直线的斜率为tan（匕）、、Xj-X3
点（X],*）与点（X3,匕）两点连线的长度：
如=/（匕一匕）2+（乂|一乂3）2;
所以这两点在一个角度为。
上的探测器上的的长度：
L!
3（6＞）'=（cos（arctan（71"^）-6＞））xJ（匕一匕―+（X】一乂3淀；乂1-乂3、、
椭圆中心点在探测器上的投影长度可以通过附件2给出的数据来拟合出一个函数来表示它：
W）=K'；
根据长度关系可得：
匕'=匕'+顷）';
即:
因为洌度并不是连续的，而是180个角度数据那么可化为：
F3（Q）,[•=1,2,3,179,180
综上所述，可以建立旋转中心模型，模型如下所示：
180180
M心Z（F3（Q）/180）-「3（%）2
j=li=l
s.t.{F3S=F\（Q）-
用MATALAB求解得：
X3=40.7543；
匕=56.2520；
E（6>,X71.1125,1=1,2,3,.180
即旋转中心坐标为（40.7543,56.2520）（相对于正方形左下角建立的直角坐标系为参考系的）。
5.2对于问题二的求解
5.2.1问题二的求解步骤
（1）求出投影坐标与原图像坐标之间的一个转化关系；
（2）假设x射线旋转角连续时，利用滤波反投影法求出经过滤波后的射线投影在。
〜180度范围的衰减系数的叠加，求出衰减系数分布函数〃（x,y）；
（3）将衰减系数分布函数#（x,y）进行离散化，利用问题一求解出的180个角度进行累加，得出待重建图像的像素点，重建出CT图像；
5.2.2问题二模型的建立与求解
5.2.2.1投影坐标与原图像坐标之间的转化
图七投影坐标与原图像坐标之间的转化关系图
由投影与几何知识可知：
（1）
（2）
p=xcosO+ysin。
；
q=-xcosO+ysin。
；
其中。
为射线的方向（x射线与原图像纵坐标的正方向的夹角）;
5.2.2.2基于滤波反投影图像的重建模型
（-）假设］射线方向角度的变化是连续的
以CT系统旋转中心为原点，X射线方向为纵坐标（g轴），垂直于X射线方向的为横坐标（〃轴）建立直角坐标系，画出I射线平行投影曲线，如图
假设七是CT系统扫描过程中得到的投影数据，5,可知乙在数值上等于1射线方向上衰减系数的线积分，由于该积分与X射线穿过模板的厚度有关，对于不规那么的图形来说，不同的投影方向或X射线离旋转中心的距离不同，尤射线穿过
模板的厚度那么不同，所以乙的大小与X射线的角度。
以及投影射线到旋转中心的距离P有关。
设f（p,q）是基于poq坐标（对应旋转后的坐标系）时，（P，q）点处的衰减系数，#（x,y）基于xoy坐标（对应原图形的坐标系）时，（x,y）点处的衰减系数，那么有：
L（p,O）=Lf（p,q）dq；
对L（p,e）进行一维傅里叶变换得：
0）=「L（p,0）e~j2mvpdp；
然后（3）代入到（4）得：
L（w，O）=LLf（P，q）e「fdpdq；
结合
（1）
（2）可得
l（w,o）=n:
"（尤,必）/5心COS"W）d裁y；
再将jn（x.y）进行二维傅里叶变换，可得：
F（s,t）=仁引埼心,y）厂j2兀（sx+ty）dxdy；
其中：
\s=wcos3
|/=wsinQ
将（8）代入（7）中可得：
F（wcos6>,wsin6>）=L（w,3）；
对（7）进行傅里叶逆变换可得：
心y）=匚匚F（s,t）『5）dsdt；
（3）
（4）
（5）
（6）
（7）
（8）
（9）
（10）
结合（8）得:
ds/dwds/dOdt/dwdt/dO
=了匚％3,。
）。
眼心cosm，s"腥湖9；
由对称性又可知:
L（p,6+兀）=L（一p,0）；
L（w,6+几）=乙（电9）；
有此可得:
"（x,y）=「「以也。
冲斯印我+网展加如心；
该函数为图像的衰减系数分布函数，表示在（x,y）处，所有经过滤波后的射线投
影在0〜180度范围的衰减系数的叠加。
（二）将衰减系数分布函数//（x,y）进行离散化
当假设X射线
升级会员 