最终版本 2017高教社杯全国大学生数学建模竞赛优秀范文

CT系统参数标定及成像问题研究

摘要

CT机扫描部分主要由X线管和不同数目的控测器组成，用来收集信息。X线束对所选择的面层进行扫描，其强度因和不同密度的组织相互作用而产生相应的吸收和衰减。

[1] 探测器将收集到的信息经过一系列的转变，最后经过计算机的储存和处理，得到CT值可以排列成数字矩阵。

通过对题目所提供材料进行分析，提出了较为合理的假设，对各组附件数据进行了拟合处理制成各种图像并分析说明，且建立模型来求解CT系统拟合处理问题。

在对问题一的分析中，对附件一模拟实体立体化建立模型Ⅰ，并对数据进行处理及排差，假设载物台在理想状态下是水平并与探测器无偏差，而且不考虑机械系数或各种问题的情况下，建立起了一个模拟CT系统的仪器。运用数学几何知识作图，通过建立相似图形（模拟CT系统运行）等比例来确定几个系统参数之间的关系（CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向）。在对建立的模型Ⅰ进行改进的基础上，对附件2进行拟合处理建立模型Ⅱ，利用数学中的傅里叶变换算法等比对图2模板示意图进行平面配对。借助数学算法和MATLAB软件，对附件中所提供的数据进行了筛选，去除异常数据，对残缺数据进行适当补充，并随机抽取了其中几组数据对理论结果进行了数据模拟，结果显示，理论结果与数据模拟结果吻合。

在对问题二的分析中，对附件3模拟建立模型Ⅲ。利用上述CT系统得到的某未知介质的接受信息还有结合问题一所得到的标定参数，通过建立相似图形等比例来确定几个系统参数之间的关系（CT系统旋转中心在正方形托盘中的位置、几何图形以及该吸收率等信息）。借助数学算法和MATLAB软件，利用图3所给的10个位置，对附件4中所提供的数据（对附件4模拟建立模型Ⅳ）进行了筛选，去除异常数据，对残缺数据进行适当补充，并随机抽取了其中几组数据对理论结果进行了数据模拟推测其的吸收率。

在对问题三的分析中，对附件5模拟建立模型Ⅴ。利用上述CT系统得到的某未知介质的接受信息还有结合问题一所得到的标定参数，通过建立相似图形等比例来确定几个系统参数之间的关系（CT系统旋转中心在正方形托盘中的位置、几何图形以及该吸收率等信息）。借助数学算法和MATLAB软件，利用图3所给的10个位置，进行了数据模拟推测其的吸收率。

在对问题四的分析中，借助数学算法和MATLAB软件，分析问题一中参数标定的精度和稳定性，并借助问题一的条件设计出新的模板、建立所对应的标定模型，以改进精度和稳定性。

关键词：数字矩阵拟合处理傅里叶变换算法平面配对标定参数吸收率

1.问题重述

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）

图1. CT系统示意图图2. 模板示意图（单位：mm）图3. 10个位置示意图

2．问题分析

问题一：

本题的误差主要在于工业CT系统的安装过程。CT系统安装时往往存在误差，从而影响成像质量，而且探测器自身的缺陷，载物台自身的不足等，通通都会成为CT系统测量的不完全的因素，所以在分析问题时，模型是假定处于理想状态下去根据数据以及利用几何知识来实现CT系统参数标定。

问题二：

根据分析，载物台很有可能会是倾斜的，而载物台可以等效为探测器投影绕一个点在竖直方向旋转。由附录一可得知其图像，又由其图像可以推测其误差来源，从而可以进行修正，最终可以给出准确的CT系统参数标定。

问题三：

根据分析，载物台是倾斜的，探测器也可能本身有缺陷，此为机械误差，而载物台可以等效为探测器投影绕一个点在竖直方向旋转。由附录可得知其图像，又由其图像可以推测其误差来源，从而可以进行修正，最终可以给出准确的CT系统参数标定。

3．模型假设

（1）假设数据来源真实有效，能够反映实际情况。

（2）各组数据在相同条件下进行，无其他人为因素影响。

（3）题中有并未明确给出但又容易得到的数据，若可测量，不列入考虑范围。（4）在数据的采集过程当中，我们认为CT系统的自身参数不会变化。

（5）在此实验中光的衍射现象不明显，在假定光以直线传播的情况下进行数据处理。

4．符号说明

5．建模过程

假设载物台在理想状态下是水平并与探测器无偏差，而且不考虑机械系数或各种问题的情况下，利用PROE建立起了一个模拟CT系统的仪器。

下图为利用附件一的数据以及MATLAB做成的附件1散点图像。

图2

对比分析：其与图2的形状大致吻合

说明附件一数据基本无偏差，完成排差工作。

2、对附件2模拟实体立体化建立模型Ⅱ

下图为利用附件2的数据以及MATLAB做成的附件2立体图像。

模型求解

（1）运用数学几何知识作图，通过建立相似图形（模拟CT系统运行）等比例来确定几个系统参数之间的关系。

①CT系统旋转中心在正方形托盘中的位置

[1]平行投影（傅里叶变换算法）【2】

将模板看成图像函数f(x,y),穿过f(x,y)的一条线称为射线，f(x,y) 沿某一射线的积分称为射线积分,而射线积分的集合则组成投影。若从坐标原点向射线作一垂线,以此垂线作为新坐标的一个轴t,并构成新坐标系(t, s)坐标仅是(x,y)坐标系旋转θ角的结果,二者存在下列变换关系:

因此射线积分可表达为：

其中射线方程可写为

其中t是原点到该射线的垂直距离。

傅里叶变换进行二维图像的重建

设错误！未找到引用源。为一二维图像，其傅里叶变换为

设该图像在x轴的投影为：

现求其傅里叶变换

通过傅里叶变换，此公式为傅里叶变换的中心剖面。

当满足时，只需对投影的傅里叶变换进行反变换，就可以得到重建的模型。

②计算探测单元之间的距离

当512个探测器处于模板的正上方时，可以利用公式

S：探测单元之间的距离 L：模板收到射线的相对长度

n：接收到信息的探测器的个数

从模板示意图可以分析出，当512个探测器处于模板正上方时，某一个探测器必定有最大吸收率，通过excel的新建规则，可以找出最大的吸收率，位置是EU223，数值为141.7794。

即EU这一列的数据为我们假设模板处于正上方时512个探测器的吸收信息，通过统计可知，模板椭圆的短轴对应的探测器个数为29，模板圆的直径对应的探测器的个数为109。

(1) (2)

(3)

联立上述3式得s=0.2755mm。

①CT系统使用的X射线的180个方向

通过MATLAB 对数据进行模拟图像

（2）在对建立的模型Ⅰ进行改进的基础上，利用对称投影相关法还有数学中的傅里叶变换算法等比对图2模板示意图进行平面配对。

5.2.问题二的建模与求解

模型建立：

1、对附件3模拟建立模型Ⅲ

下图为利用MATLAB做成的附件3图像

2、对附件4模拟建立模型Ⅳ

下图为利用MATLAB做成的附件4散点图像

结合图三，进行对比

发现图形与数据吻合，通过MATLAB,数学知识及图形比对大致可以得出其吸收率。

5.3.问题三的建模与求解

模型建立

1、对附件5模拟建立模型Ⅴ

下图为利用MATLAB做成的附件5图像

模型求解

利用MATLAB做成的附件4散点图像图三

4.1.问题四的建模与求解

模型建立

1、CT成像原理

Radon变换将图像投影到某个角度的放射线上。

2、公式表达【3】

Radon逆变换的公式为：

投影射线方程可表示为

3、模型假设

4、验证精度与稳定性

利用MATLAB对数据进行重新调整，测试其精度与稳定性，以下图片为其论证

6.模型的综合评价

模型的优点：Radon变换可用于直线检测，它的积分运算抵消了噪音所产生的误差。Radon变换可以将低质的不清楚的盲点图像恢复成原始图像。

Radon可以灵活选取不同积分路径，由此算出一系列相似变换。

模型的缺点：Randon变换的离散化是一个比较复杂的问题，在众多的离散化算法中，有些存在大量的冗余，有些虽然克服了大的冗余度，但是得到其所对应的逆变换又比较困难。其中有限Radon变换FRAT是其中比较好的离散化算法之一。有限Radon变换是有限大小的二维离散图像实现Radon变换的离散化方法。

模型的改进：为了获得更好的能量集中性，有限Radon变换（FRAT）和反变换FBP要求变换的图像均值为零，对于均值不为零的图像可以在变换前先减去均值，以保证变换前的图像均值为零；反变换回来后再加上图像均值即可恢复原图像。

模型的推广：此模型不仅可以运用在CT领域上，还可以运用在炮集记录

7.参考文献

[1] 叶宁军，《CT系统简介》，科学出版社，2009年

[2] 朱翚、王富东，《利用matlab实现二维图像傅里叶变换法》，科学出版社，2012

年

[3] 郭立倩，《CT系统的标定与有限角度CT重建方法的研究》，浙江科学技术出版社，

2003年

[4] 巩向博，《高精度Radon变换及其应用研究》，高等教育出版社，2005年

[5] 《Matlab R2016a完全自学一本通》，中国工信出版社，2006年

1.附录

附录1

function createfigure(X1, Y1)

%CREATEFIGURE(X1, Y1)

% X1: x êy?Yμ?ê?á?

% Y1: y êy?Yμ?ê?á?

% óé MATLAB óú 15-Sep-2017 21:38:09 ×??ˉéú3é

% ′′?¨ figure

figure1 = figure;

% ′′?¨ axes

axes1 = axes('Parent',figure1,...

'Position',[0.118609680865194 0.115263157894737 0.781900103723995

0.809736842105265]);

hold(axes1,'on');

% ′′?¨ plot

plot(X1,Y1,'DisplayName','?üê??êé￠μ?í?','Marker','.','LineStyle','none',... 'Color',[0 0.447 0.741]);

% ′′?¨ xlabel

xlabel('nz = 12568');

% è???ò???DDμ?×￠êíò?±￡á?×?±ê?áμ? X ·??§

% xlim(axes1,[0 256]);

% è???ò???DDμ?×￠êíò?±￡á?×?±ê?áμ? Y ·??§

% ylim(axes1,[0 256]);

view(axes1,[-0.800000000000026 90]);

box(axes1,'on');

axis(axes1,'ij');

% éè????óà×?±ê?áê?D?

set(axes1,'GridLineStyle','none','PlotBoxAspectRatio',[257 257 1],...

'XTickLabel',{'0','20','40','60','80','100'},'YTick',[0 50 100 150 200 250],...

'YTickLabel',{'0','20','40','60','80','100'});

% ′′?¨ legend

legend1 = legend(axes1,'show');

set(legend1,...

'Position',[0.619871268225592 0.856578948204979 0.116374869056882

0.0350877184616892]);

附录2

function createfigure(zdata1)

%CREATEFIGURE(ZDATA1)

% ZDATA1: surface zdata

% óé MATLAB óú 15-Sep-2017 21:46:22 ×??ˉéú3é

% ′′?¨ figure

figure1 = figure;

% ′′?¨ axes

axes1 = axes('Parent',figure1);

hold(axes1,'on');

% ′′?¨ mesh

mesh(zdata1,'Parent',axes1);

view(axes1,[-37.5 30]);

grid(axes1,'on');

附录三

function createfigure(zdata1)

%CREATEFIGURE(ZDATA1)

% ZDATA1: surface zdata

% óé MATLAB óú 15-Sep-2017 21:48:07 ×??ˉéú3é

% ′′?¨ figure

figure1 = figure;

% ′′?¨ axes

axes1 = axes('Parent',figure1);

hold(axes1,'on');

% ′′?¨ mesh

mesh(zdata1,'Parent',axes1);

view(axes1,[-37.5 30]);

grid(axes1,'on');

附录四

function createfigure(X1, Y1, S1, C1)

%CREATEFIGURE(X1, Y1, S1, C1)

% X1: scatter x

% Y1: scatter y

% S1: scatter s

% C1: scatter c

% óé MATLAB óú 15-Sep-2017 21:49:39 ×??ˉéú3é

% ′′?¨ figure

figure1 = figure;

% ′′?¨ axes

axes1 = axes('Parent',figure1);

hold(axes1,'on');

% ′′?¨ scatter

scatter(X1,Y1,S1,C1);

附录五

function createfigure(zdata1)

%CREATEFIGURE(ZDATA1)

% ZDATA1: contour z

% óé MATLAB óú 15-Sep-2017 21:51:38 ×??ˉéú3é

% ′′?¨ figure

figure1 = figure;

% ′′?¨ axes

axes1 = axes('Parent',figure1);

hold(axes1,'on');

% ′′?¨ contour

contour(zdata1);

box(axes1,'on');

axis(axes1,'tight');

% éè????óà×?±ê?áê?D?

set(axes1,'BoxStyle','full','Layer','top');