CT 系统参数标定及成像
摘要: 本文运用 MATLAB 等工具对已给出的数据进行分析和处理,通过反
射投影算法,等比例转换法, radon 变换和 iradon 变换,还原 180 次扫描信息和
图形信息。
对于问题 1,通过 radon 变换法,在 MATLAB 中得出该介质以正方形托盘
左上角为原点的坐标系下的位置分布图, 然后根据题目中已经给出的介质物体实
际图,以椭圆圆心为原点建立直角坐标系, 得出两个坐标系之间的比例关系, 通
过位置与长度的等比例变换得出旋转中心在正方形托盘中的坐标为(
-8.7755,
6.1697),通过观察附件 2 发现存在探测器接收到的非零信号个数稳定在 28 个,
对比小圆的直径得出探测器的间距为 0.2857mm,探测器接收到的非零信号个数
与角度曲线没有发生突变,且最高点与最低点横坐标相差
90 次,可以认为每次
旋转 1 度,初始位置与坐标系 X 轴正方向夹角为 29 度。
对于问题 2,通过使用 iradon 变换,得出了投影重建结构的解,对附件 3 中
某未知介质的投影数据进行滤波反投影重建运算, 实现从其它空间向图像空间进
行转换的过程,最终通过 MATLAB 运行结果获得该未知介质模型的重建图像,
得出该未知介质在正方形托盘中的几何形状和位置信息, 然后采用比例变换的方
式,根据 10 个点的位置和相对于实物图位置, 得出这 10 个位置介质点的吸收率
结果。
对于问题 3,采用与问题 2 相似的方式,利用 MATLAB 中的 iradon 算法,
根据附件 5 中提供的另一未知介质的吸收信息, 通过反投影重建可以得到该未知
介质的位置, 形状和吸收率等信息, 同样采用等比例变换的方式, 根据点的位置
和相对于实物图的位置,得出这 10 个位置点的吸收率结果。
对于问题 4,通过对已经给定的数据进行分析,用 iradon 验证扫描次数对成
像质量的影响,在不同滤波环境下比较成像质量,分别对 18,36,90,180个角度投
影进行观察和分析,能够得出随着投影角度个数的增加,图像的重影越来越少,
也即是稳定性和精确度越来越高。运用 shepp-lagon模型重新优化模型。
关键词: 反射投影重建; MATLAB 软件; radon 变换; iradon 变换;比例变
换;成像质量;
1.问题的重述 :
CT 是一种利用样品对射线能量的吸收特性对生物组织和工程材料的样品进
行进行断层成像的技术, 由此可以获取样品的内部结构信息。 本次使用的二维 CT
系统中平行入射的 X 射线垂直于探测器的平面,并且探测器单元可视为等距离
排列的接收点, X 射线的发射源和接收源的相对位置不变, 整个系统绕着某个固
定的旋转中心逆时针旋转 180 次,每个 X 射线方向,具有的 512 个等距离单位
元探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量,
并经
过增益等处理后得到 180 组接收信息。 但是 CT 系统在安装时往往存在一定的误
差而给测量质量带来一定的影响, 本次借助于已知结构的样品标定 CT 系统的参
数,题目要求建立相应的数学模型和算法解决以下四个问题:
(1)根据正方形托盘上的两个均匀介质模板的几何信息以及其所反应的吸
收率和接收信息,确定 CT 系统旋转中心在正方形托盘中的具体位置,探测器单
元之间的距离以及该 CT 系统所使用的 X 射线的 180 个方向;
(2)利用 CT 系统得到的某未知介质的接收信息以及 1 中的标定参数,确
定未知介质在正方形托盘中的位置,几何形状和接收率,以及图
3 给定的 10 个
位置处的吸收率;
(3)利用 CT 系统得到的另一未知介质的接收信息以及 1 中得到的标定参
数,确定 该未知介质的信息,以及图 3 给定的 10 个未知处的吸收率;
(4)分析 1 中参数标定的精度和稳定性,自行设计新模板、建立对应的标
定模型,改进标定精度和稳定性,给出理由;
2.问题的分析:
CT 技术是一种依据外部投影数据重建被探测物体的内部结构的无损检测技
术,具有高精度,高效率,无损检测等优点,近年来被广泛应用于生物医学,航
空航天,材料加工,组织探伤等领域,图像重建算法是 CT 技术的核心之一,重
建算法的优劣直接影响了 CT 检测性能的好坏,而对于未进行图像校准的 CT 检
测系统来说, 确定其旋转中心和每次旋转的角度又极为重要, 下边的几个问题就
是解决 CT 旋转系统的中心,每次旋转的角度,通过建立 CT 成像算法,来进行
模型检验和优化求解问题。
(1) 对于问题 1,通过对附件中所给出的相应数据进行分析, 附件 1 中是
将图 2 模板示意图利用建立二维坐标系,将原本连续的模板划分成独立的小单
元,每个单元都能吸收 X 射线,我们可以把每个单元看成一个个像素, X 射线从
一个方向射向模板,,穿过模板后被探测器检测到,而每个探测器接收到的信号
强弱各有不同,因此可以简单分析认为,探测器的信号和被接受的
X 射线经过
模板的长度有关, CT 围绕旋转中心旋转进行探测,每次旋转过一定的角度,而
模板并不是一个完全对称的物体, 在各个方向上的长度各有不同, X 射线通过的
长度不同, 2 中已经显示每次旋转得到的探测其能量和旋转角度关系,通过分
析某一固定位置的探测器接收到不同能量的变换来求得每次旋转的角度, 探测器
的间隔可以根据模板的长度和探测到 X 射线能量变化的接收器个数来确定。
(2)对于问题 2,需要使用在问题 1 中所求出的 CT 系统旋转中心在正方形
托盘中的位置,以及探测器单元之间的距离和 CT 系统中 X 射线的方向这些结
果,附件 3 给定了未知介质的接收信息数据, 通过对其中的数据进行分析, 使用
与问题 1 互逆的方式, 通过使用算法求出正方形托盘上未知介质的相对位置, 几
何形状和吸收率, 另外还给出了正方形托盘中的 10 个点的位置, 来得出这 10 个
点的吸收率,在解决问题 1 时,是使用吸收率得到相对位置,因而求 10 个点的
吸收率也是使用进行问题 1 的逆向求解过程。
(3)对于问题 3,可以采用与问题 2 相似的方法,也是对问题 1 的过程进
行逆向求解,从而得出该未知介质的相关信息,以及另外
10 个位置点的吸收率
的大小。
(4)对于问题 4,通过结合前面 3 个问题的结论与实现方法,自行设计一
种模板和相应的模型, 对于处于不同位置的介质点, 其吸收率也会不尽相同, 可
以仿照附件 1 中的模板介质的吸收率和附件 2 与附件 3 的介质相关信息, 先设定
一组吸收率或者介质的相关信息, 使用问题 1 的逆向求解过程, 得出模板的位置
和几何形状信息,并求出问题 1 和问题 4 中相应的参数标定误差和稳定性的大
小。
3.模型的假设和符号的说明
3.1 模型的假设 :
(1)假设 X 射线穿过介质时, 由于 X 射线与介质的原子相互作用过程中的
散射,衍射等对 X 射线强度的衰减的影响可以忽略不计,也即是 X 射线衰减只
与介质吸收有关;
(2)假设整个 CT 发射 -接收系统绕着旋转中心旋转的过程是匀速进行的;
(3)假设入射的 X 射线的能量分布是均匀的,且每次发射源发射的 X 射线都
是相同的;
(4)假设介质是等密度分布的,单位质量介质中的电子数目相同;
3.2 符号的说明
I
I0
P
μ
Δl
f(x,y)
n
pi
符号说明表
dij
λ
入射光强
出射光强
点 i→j 的距离
两 直 角 坐 标 系
间的比例系数
探测器探测到的
h(x,y)
待 变 换 目 标 函
能量
衰减系数
单位方格宽度
衰减系数函数
扫过圆的探测器
个数
第 i 个点的坐标
q
w
T
β
数
傅 立 叶 变 换 原
函数
吸收强度
灰度值
w 与 t 间的比值
4.模型的建立与求解:
4.1 问题 1 的模型建立与求解:
4.1.1:求初始位置角度和探测器单元间距离:
由朗伯比尔定律可知:用一束强度为 I0 的 x 光照射物体,如果物体内部均匀,
已知物体对于 x 光的衰减系数为 μ,那么出射光强 I 与入射光强 I0 存在如下关
系:
探测器探测到的能量:
I = I 0?- μc
p = In(I 0 ∕I) ,可知, p = - μc
在沿某个方向入射的 X 射线,探测器探测的能量与 X 射线经过的物体的长度
成正比。
f i 为 x-y 平面的衰减系数函数,然后将物体划分成大小为 ?l ×?l 均匀的方格,
使其离散化。
I0
f2
f3
f1
?l
fn
I
图 4-1 朗伯比尔定律示意图
探测器接收到的能量
p = ∑?????l=ln
I 0
I ,
当?l 趋向于无穷小时,得到线积分
p = ∫????????
针对假设,我们得知在某个方向上, p 可以看做在该方向上所有射线的投影集
合,通过求解一定角度的 p 来求解物体内部平面各点的衰减系数函数 f(x,y) 。
图 4-2-1 模板示意图(单位: mm)
图 4-2-2 反射投影模板结构图·
利用 MATLAB 绘制的附件中模板图中椭圆和圆的方程:
椭圆方程为:
圆方程为:
x2
225
+
??2
1600
= 1
(x - 45)
2
+ ??2 = 42
选取椭圆短轴所在直线为 x 轴,长轴所在直线方向为 y 轴
平行束反投影重建算法 (Filter back projeCTion,abbr. FBP)
傅立叶中心切片定理:
假设 f(x,y) 是待重建物体的密度函数, p(ti, α)为 f(x,y) 在角度 α时的平行束投影,
有傅立叶中心切片定理的数学表达式为:
F??= ??(??,??)
F1 表示一维傅立叶变换 F( ρ, α)是二维傅立叶变换的极坐标表示。
Radon变换就是将数字图像矩阵在某以指定角度射线方向上做投影变换, 例如,
二维函数的投影就是在其指定方向上的线积分, 在垂直方向上的线积分就是在 x
轴上的投影;在水平方向上的二维线积分就是在 y 轴上的投影。
设直角坐标系( x,y)转动 α角后得到旋转坐标系( x′,y′),由此得知:
x′=xcos α+ysin α
p(x ′,??)为原函数 f(x′,y ′) 的投影 f(x,y) 沿着旋转坐标系中 x′轴 α方向的线积分,
根据定义公式知其表达式为:
p(x ′,??)=?
∞
∞
f(x,y)
??( ??????????
+ ??????????
- ??)????????, 其中 0 ≤??≤ π
我们构建接收器接收到的信号和经过 Radon变换得到的图形投影变换模型:
接收到衰减信号的接收器个数 t(ti)正比于穿过图形的个数,如图所示:
P=-μ*p(x′, ??)
图 4-3 接收到衰减信号的接收器个数
从 A 题附件得知在 DE 列之后,图像分为两部分;分析可知,对于模型有一
个椭圆和一个小圆,根据提出的模型,接收器接收到的信号
p>0 的区域也分为
两个部分,因此,假设两者之间存在对应关系,选取从
A 列至 FX 列、从 111
行数据利用 matlab 统计存在 p>0 的个数,如图所示:
图 4-4
111 行能接受到非零信号的探测器个数与转动次数示意图
可知从 DE 列到最后 72 列探测器信号 p>0 的个数在一定范围内趋于一致,利
用 MATLAB 求其平均值, n 为个数, n?=∑????(??> 0)
求得 n?=28.8333, 取整n? = 29,对应于小圆宽度为 2*r
d=
2???
n?-1 =0.2857mm
对 A 题附件 2 进行 n(p>0)统计,利用 MATLAB 得出如下图形:
图 4-5 能接受到非零信号的探测器个数与转动次数图
通过图形能够得出最大值为 289,对应的转动次数分别为 58、59、60、61、
62、63、64、64,取其均值作为存在最大值的转动次数,均值为 61,最小值为
137,对应的转动次数为 150、152,取其均值作为存在最小值的转动次数,均
值为 151,根据图像是平滑的曲线,不存在任何突变,且不存在两组完全相同
的探测器数据,因此可以假设 CT 为均匀旋转,且每次旋转 1 度,则初始位置
为 29 度。
4.1.2 CT 系统旋转中心位置:
判断本 CT 系统的旋转中心, 我们通过 radon 变换法,在 MATLAB 中模拟出介
质物体的中心点的位置如下图所示,并且能够得出相应点的坐标位置: