2017年国赛高教杯数学建模A题CT系统参数标定及成像解题全过程文档及程序

2017年国赛高教杯数学建模

A题 CT系统参数标定及成像

  CT(Computed Tomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。一种典型的二维CT系统如图1所示，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。
  CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。
  请建立相应的数学模型和算法，解决以下问题：
  (1) 在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，相应的数据文件见附件1，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对应于该模板的接收信息见附件2。请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。
  (2) 附件3是利用上述CT系统得到的某未知介质的接收信息。利用(1)中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率，相应的数据文件见附件4。
  (3) 附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。
  (4) 分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。
  (1)-(4)中的所有数值结果均保留4位小数。同时提供(2)和(3)重建得到的介质吸收率的数据文件（大小为256×256，格式同附件1，文件名分别为problem2.xls和problem3.xls）

整体求解过程概述(摘要)

  本文针对CT系统参数标定及成像问题，建立了基于最小二乘法的CT系统标定模型，求解了CT系统的旋转中心坐标、探测器单元间距以及X射线的180个方向；建立了基于傅立叶变换的图像重建模型，求解了未知介质的吸收率分布；建立了基于滤波反投影法的图像重建模型，求解了噪声较大的未知介质的吸收率分布；设计了新型标定模板并建立了对应的标定模型，解决了改进标定精度和稳定性的问题。
  针对问题一，建立了基于最小二乘法的CT系统标定模型，求解了CT系统的标定参数。首先，通过最小二乘法，将180组接收信息拟合成正弦函数，根据几何关系联立方程组求解出180个旋转角度，以及探测器单元间距𝑑=0.2758。其次，通过最小二乘优化、变步长遍历搜索的方法，求得旋转中心为(-9.2513,5.6783)。最后，基于所得标定参数及已知模板的吸收率分布，求解标定模板的接收信息，进行残差分析，检验了标定模型的合理性。
  针对问题二，建立了基于傅立叶变换的图像重建模型，求解了未知介质的吸收率分布。首先，依据 Radon 变换及中心切片定理，对投影数据进行一维傅立叶变换，得到定义在傅立叶空间的极坐标网格。其次，在频域上对投影进行积分，再次，直接通过二维傅立叶反变换得到重建图像，其中10处吸收率分别为0.0000、0.0000、1.0060、0.0000、1.0081、0.0081、0.0000、0.1008、0.0081、0.0081。最后，对标定参数进行了灵敏度分析。
  针对问题三，建立了基于滤波反投影法的图像重建模型，求解了噪声较大的未知介质的吸收率分布。首先，采用问题二的方法，发现重建图像噪声较大，改用滤波反投影算法去除重建图像的伪影。再次，选择中值滤波器进行滤波卷积，采用线性内插的方法进行插值，反投影累加后得到重建图像，其中 10处吸收率分别为 2.0323、1.2863、2.0081、1.6370、2.0484、2.0081、0.0000、2.0343、1.5363、0.0403。最后，对标定参数进行了灵敏度分析。
  针对问题四，设计了新型标定模板并建立了对应的标定模型，解决了改进标定精度和稳定性的问题。首先，对原有模板进行误差分析，得出对椭圆边缘处接收信息的插值拟合会导致旋转角度出现误差的结论。其次，为避免类似误差，设计圆形组合模板，并建立对应标定模型。再次，对比基于新旧模板的不同标定参数的图像重建模型，求得新旧模板的峰值信噪比分别为50.3614、52.2537。最后，求解基于不同旋转角度个数的归一化均方距离，新模板重建图像精度趋于稳定的旋转次数小于旧模板，进而验证了基于新模板的标定模型精度更高，更稳定。

模型假设：

  1. 假设射线源中点与转台旋转中心连线垂直于探测器阵列；
  2. 假设射线源的衰减仅与被介质吸收有关，并未产生散射；
  3. 假设增益处理是线性增大的。

问题分析：

  问题整体思路较为清晰，层层递进。问题一是根据已知标定模板的吸收率及其接收信息对CT系统进行正向标定的问题。问题二、三建立在问题一求解的标定系数基础之上，是依据不同的接受信息对未知介质的吸收率分布进行反向重建的问题。分析问题二、三的接受信息可发现问题三图像噪声偏大，故对问题三的求解大致思路与问题二类似，但需进行进一步的降噪处理。问题四需对原有标定模型进行精度和稳定性分析，并设计新的模板及标定模型，进而提高精度及稳定性，属于优化问题。
  问题一的分析
  首先，分析题目条件可知，均匀模板的吸收率分布及对应的接受信息已知，且模板质地均匀，故射线衰减公式中的对吸收率分布进行的线积分可简单变为介质厚度，故可对接收信息的表达式进行求解。 其次，要将离散的接受信息连续化，可考虑通过最小二乘法对其进行拟合， 初步分析为两个近似正弦曲线，分别对应图中的椭圆和圆，其中正弦曲线与X轴相交的情况反应了平行光束与椭圆相切的两条切线距离的大小。由于两切线之间的距离与旋转中心[1]、旋转角度存在一定的几何关系，列出关系式，联立方程可解出各个旋转角度的值。 再次，基于上述拟合曲线及几何关系式，可求得多个探测器单元间距及旋转中心的坐标。针对探测器单元间距的求解，出于对使用信息完整性的考虑，可通过多次求解方程组，并对多个解取平均得到最终答案。针对旋转中心的求解，由于探测器中垂线衡过旋转中心，且旋转角度已求解得到，故可列出不同角度下中垂线的直线方程，两两联立即可求出旋转中心。出于保留信息完整性，同时避免算法过于复杂的考虑，故均匀抽取30组直线，得到元素个数为345的点集。可通过最小二乘优化[2]，建立各点到中心点距离平方和最小的目标函数，可通过变步长遍历搜索的方法，求解旋转中心。 最后，由于接受信息的表达式可被求解，且介质吸收率分布已知，故可对已知模板的接受信息进行模拟，可与题中所给条件进行比对，进行残差分析，即可检验模型的合理性。
  问题二的分析
  首先，分析题目条件可知，由于在第一问已对CT系统进行了标定，故可根据题中所给的接收信息，对图像进行重建。对图像重建的常用方法有傅里叶变换法、滤波反投影法，由于部分接收信息的拟合后的图像较为平滑，故采用算法较为简单的傅里叶变换法。 其次，在推导非均匀介质中衰减公式时，可对非均匀物体进行微元化考虑，进而对射线贯穿部分的吸收率分布进行线积分。可根据Radon变换[3]及中心切片定理[4]，得到多个角度下的投影后，通过投影反变换可求出被积函数即线衰减系数函数，从而得到对应的线衰减系数分布的图像，重建出问题二中非均匀物体的形状及在正方形托盘上的位置。 再次，为得到题中所要求的10个点的吸收率，需要得到吸收率与图像灰度值之间的关系。由于问题一中模板条件已知，运用问题二中的算法可重建出问题一中的模板图像，得到已知模板图像与重建模板图像之间的灰度值关系后再据此关系得问题三图像的10个吸收率。 最后，可针对三个不同标定参数，进行固定比率的修改，并计算参数修改后重建图像吸收率分布的变化率，可据此对各标定参数进行灵敏度分析。
  问题三的分析
  首先，问题三与问题二提供的条件并无二致，然而通过观察附件中的接收信息值，可发现数据在表格中的分布更加复杂，很可能导致运用问题二中的傅