2017年大学生数学建模A题CT系统标定成像论文

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 入射光强 dij 点 i j 的距离 I0 出射光强 两直角坐标系 间的比例系数 P 探测器探测到的 能量 h x y 待变换目标函 数 衰减系数 q 傅立叶变换原 函数 l 单位方格宽度 w 吸收强度 f x y 衰减系数函数 T 灰度值 n 扫过圆的探测器 个数 w 与 t 间的比值 pi 第 i 个点的坐标 I0 f1 f2 f3 fn l 4 模型的建立与求解模型的建立与求解 4 1 问题问题 1 的模型建立与求解 的模型建立与求解 4 1 1 求初始位置角度和探测器单元间距离 求初始位置角度和探测器单元间距离 由朗伯比尔定律可知 用一束强度为 I0的 x 光照射物体 如果物体内部均匀 已知物体对于 x 光的衰减系数为 那么出射光强 I 与入射光强 I0存在如下关 系 I I0 c 探测器探测到的能量 p In I0 I 可知 p c 在沿某个方向入射的 X 射线 探测器探测的能量与 X 射线经过的物体的长度 成正比 fi为 x y 平面的衰减系数函数 然后将物体划分成大小为 l l均匀的方格 使其离散化 图 4 1 朗伯比尔定律示意图 探测器接收到的能量 p l lnI0 I 当 l趋向于无穷小时 得到线积分 p 针对假设 我们得知在某个方向上 p 可以看做在该方向上所有射线的投影集 合 通过求解一定角度的 p 来求解物体内部平面各点的衰减系数函数 f x y I 图 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 变换得到的图形投影变换模型 P p x 接收到衰减信号的接收器个数 t ti 正比于穿过图形的个数 如图所示 图 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 中模拟出介 质物体的中心点的位置如下图所示 并且能够得出相应点的坐标位置 图 4 6 介质物体位置分布示意图 通过结合题目中给出的椭圆形模板和圆形模板 进行坐标系之间的等比例变 换 通过建立相应的平面直角坐标系 可以分别得到椭圆圆心 圆的圆心在该平 面直角坐标系里的坐标分别为 P1 273 6658 291 6983 和 P2 415 2547 372 0849 假设旋转中心点 P3 256 0000 256 0000 由平面内两点之间的距离公式 d Xi Xj Yi Yj 可得 d12 162 8171 d13 39 8303 d23 197 0730 由 P1和 P2点的坐标可得通过点 P1和 P2的直线的方程为 Y1 0 5677 X 136 3382 设 P3点与直线 Y1相交于 P4 点 由平面直线之间的垂直关系可得 P4 244 9793 275 4129 可得通过P3和P4的直线方程为Y2 1 7615 X 706 9440 从而可得出交点P4 与椭圆圆心PI和旋转中心 P3之间的距离分别为d14 31 7511 d34 22 3230 在题目中所给的模板示意图中 建立以椭圆中心为原点 以椭圆圆心和圆圆心 连线为 x 轴的另一个平面直角坐标系如下 图 4 7 以椭圆圆心为原点的正方形托盘平面直角坐标系 可得在此平面内椭圆圆心 Q1 0 0 圆的圆心 Q2 45 0 二者之间的距离 D12 45 由此可得两个平面直角坐标系之间的相对比例系数 D12 d12 记在该坐标系 下的旋转中心坐标为 Q3 x0 y0 由坐标系之间的对应关系发现 Q3位于第二象限 之内 其横纵坐标分别为 x0 d14 8 7755 y0 d34 6 1697 可得该 CT 系统的旋转中心坐标为 Q3 8 7755 6 1697 以在正方形托盘平面 内 以椭圆圆心为原点建立直角坐标系 可得旋转中心点 Q3相对于椭圆和圆的 位置如下图所示 图 4 8 旋转中心的相对位置示意图 4 2 问题问题 2 的模型建立与求解的模型建立与求解 4 2 1 iradon 变换求未知介质的形状与结构 变换求未知介质的形状与结构 在图像处理中 我们一般会将 iradon 函数用到 X 射线断层扫描中 irodon 变换 是 radon 变换的逆过程 给出了投影重建结构的解 下面利用傅里叶反变换计算 iradon 函数的表达式 h x y 为需要进行变换的目标函数 用 2 D 傅里叶反变换写成极坐标形式为 dpqpjqtH 2exp q d y h x 0 其中方括号内的表达式为 q H q t 的 1 D 傅里叶反变换的表达式 结合傅里 叶变换的卷积定理 将 q 的 1 D 傅里叶反变换写成为 H 1 q H 1 q sgnq H 1 j2 q H 1 sgn q j2 利用微分性质可得 F 1 j2 p 经过整理可得 radon 反变换的表达式为 h x y 1 2 0 d htpR 1 2 1 p 根据上述的 radon 反变换的表达式 也即是 iradon 变换 通过使用 MATLAB 中的 iradon 函数实现反投影重建算法 对附件 3 中已经给出的某未知介质的投 影数据进行滤波反投影重建运算 从而从其它空间向图像空间进行转换的过程 最终通过 MATLAB 运行结果获得该未知介质模型的重建图像如下 图 4 9 未知介质模型重建图像示意图 从上图中能够直看出该未知介质在正方形托盘中的位置 也能够能够较为直观 看出该未知介质在某个空间角度下的结构图 4 2 2 求未知介质的求未知介质的相关相关信息 信息 通过使用 iradon 函数将附件 3 中的接收信息矩阵 W1 重建 得到图 然后将该图片的值矩阵导入到到 excel 图表中 基于此原理 对附件 2 采用相 同的方式 得到各自的灰度值矩阵 对于附件二 其模型介质均匀 吸收强度均 为 1 其灰度值保持在 0 5 左右 干扰值除外 即数值接近 0 的数值 根据附件 二的结果和关系 可以得到吸收强度 W 与灰度值 T 之间的相关比例为 W T 1 0 5 2 即吸收强度 W 灰度值 2 2T 因此 在去除干扰点后 把上式应用到未知介质的模型的灰度值矩阵上 可得 到对应点上的吸收率 选取部分数据如下 表 4 1 第 87 95 行 159 161 列的吸收率数据表 FC FD FE 87 0 5874 0 5893 0 5917 88 0 5879 0 5931 0 5885 89 0 5927 0 5793 0 6001 90 0 5920 0 5863 0 6027 91 0 5903 0 5917 0 5888 92 0 5865 0 5941 0 5906 93 0 5737 0 5945 0 5807 94 0 5834 0 5779 0 5820 95 0 5841 0 5795 0 5861 4 2 3 求求 10 个位置点的吸收率个位置点的吸收率 采用比例变换的方式 根据附件 4 给出的 10 个点的相对位置 图 4 10 附件 4 中 10 个点的位置分布图 相对于实物图 可在吸收强度矩阵的对应位置中找到这 10 个目标点的吸收率 如下 表 4 2 问题 2 的 10 个位置的吸收率 位置 X Y 吸收率 1 10 0000 18 0000 0 2 34 5000 25 0000 0 3 43 5000 33 0000 1 4 45 0000 73 5000 0 5 48 5000 55 5000 1 6 50 0000 75 5000 0 7 56 0000 76 5000 0 8 65 5000 37 0000 1 9 79 5000 18 0000 0 10 98 5000 43 5000 0 4 3 问题问题 3 的模型建立与求解 的模型建立与求解 4 3 1 求未知介质相关信息 求未知介质相关信息 采用与问题 2 相同的方式 利用 MATLAB 中的 iradon 算法 根据附件 5 中提 供的另一未知介质的吸收信息 通过反投影重建可以得到有滤波和无滤波条件下 该未知介质的几何形状分别为 图 4 11 有滤波和无滤波条件下未知介质几何形状图 从上面的两图中能够看出 在有滤波和无滤波条件下的未知介质的几何形状 无明显差别 4 3 2 求 10 个位置点的吸收率 采用比例变换的方式 根据附件 4 给出的 10 个点的相对位置 得出这 10 个位置 点的吸收率为 表 4 3 问题 3 的 10 个位置的吸收率 位置 X Y 吸收率 1 10 0000 18 0000 0 2 34 5000 25 0000 0 3 43 5000 33 0000 1 5470 4 45 0000 73 5000 1 4459 5 48 5000 55 5000 1 7115 6 50 0000 75 5000 3 3970 7 56 0000 76 5000 0 8 65 5000 37 0000 2 4807 9 79 5000 18 0000 0 10 98 5000 43 5000 0 4 4 问题问题 4 的模型建立与求解 的模型建立与求解 4 4 1 建立新模板与对应的标定模型 建立新模板与对应的标定模型 针对参数标定的精度和稳定性 首先对题目给出的数据进行分析 本题目是 在 180 个扫描方向上得到的扫描信息 用 iradon 验证扫描次数对成像质量的影 响 另外 根据滤波反射投影模型 在不同滤波环境下比较成像质量 进一步 分析确定参数标定的精度和稳定性 新模板采用 Shepp Logan 模型进行标定 分别得到 4 组不同个数角度投影对应的图形如下 图 4 12 不同个数角度下 iradon 变换图 通过对上面的图进行对比 18 个角度投影图中存在较多的重影 说明其稳定性 和精确度较低 36 个角度投影图中的重影比 18 个角度投影图中的重影少 90 个 角度投影图中的重影与 180 个角度投影图中的重影更少 并且通过观察可得 随 着投影角度数目的逐渐增多 所得到的投影图的重影数目越来越少 对图形稳定 性和精确性的影响也越来越小 5 模型的评价 1 将连续的二维物体碎片化 使其分散成像素点 提升了运算速度 2 模型基于在 180 个方向所得的探测器信号进行 iradon 变换 基于傅里 叶变换和反射投影算法 利用 MATLAB 处理数据 该大大提升了数据 处理能力 3 对模型中增添不同的滤波函数 提升成像的品质 将 iradon 变换得到 的灰度图类比吸收率简化算法 4 对旋转中心求解方法利用通过定点的相交直线求得 5 计算效率较低 模型处理的数据类型比较单一 参考文献参考文献 1 高飞 MATLAB 图像处理 375 例 M 北京人民邮电出版社 2015 186 190 2 马晨欣 CT 图像重建关键技术研究 D 河南 解放军信息工程大学 2011 3 张顺利 李卫斌 唐高峰 滤波反射图像重建算法研究 N 咸阳师范学院 学报 2008 8 23 4 刘明进 工业 CT 系统旋转中心定位方法研究 D 重庆 重庆大学 2014 5 毛小渊 二维 CT 图像重建算法研究 D 江西 南昌航空大 2016 6 孟凡勇 李忠传 杨民 李静海 CT 旋转中心的精确确定方法 J 中国体视 学与图像分析 2013 18 4 附录 问题 2 程序如下 1 图 4 2 2 程序 A xlsread d A 题附件 2 R F iradon A 29 208 512 imshow R 2 图 4 4 程序 clc d xlsread F 建模相关 A2 1 jishu2 zeros 180 1 jishu2 为小圆对应的个数 for i 109 180 for j 1 111 if d j i 0 jishu2 i 1 jishu2 i 1 1 end end end 111 jishu2 109 180 plot 111 jishu2 xlabel 转动次数 ylabel 非零信号探测器个数 mean 111 jishu2 求出平均值 title 111 行能接受到非零信号的探测器个数与转动次数图 ma i1 max jishu 求出最大数值和所在位置 3 图 4 5 程序 data xlsread F 建模相关 A2 1 jishu zeros 180 1 for i 1 180 for j 1 512 if data j i 0 jishu i 1 jishu i 1 1 end end end plot 512 jishu xlabel 转动次数 ylabel 非零信号探测器个数 title 能接受到非零信号的探测器个数与转动次数图 a b max 512 jishu 最大值和位置 a b min 512 jishu 最小值和位置 4 图 4 6 程序 A xlsread d A题附件 2 R F iradon A 1 512 bw im2bw R 0 1 L bwlabel bw stats regionprops L Centroid imshow L hold on for