2017数学建模A题二等奖论文

1、1 ct 系统参数标定及成像 摘要 目前，如何通过 ct 机的探测器已知的投影信息还原出物体几何形状以及内 部结构，探测已知模块得到的数据推导 ct 机的参数以及校正，是研究精密仪器 的重要方面。本文通过建立数学模型，给出了通过分析表格数据来解决上述所有 问题的方法。 对于问题一， 以正方形左下角为原点建立平面直角坐标系对方向角进行几何 抽象， 通过介质上任意一点在转动过程中在光屏上的投影点的运动轨迹都是在旋 转中心的投影点的上下波动这一几何特性，引入三角函数、 直流分量和交流分量 并与几何方法结合，建立了已知模板接收信息、吸收率、ct 参数、余弦轨迹的 傅里叶展开之间关系的综合模型。 101

2、 cos wsin w n yaanbn，tan(tan ) kk yxbag 求解模型时，将附件的探测器接收数据和带入，利用光路方向角与椭圆中心 和圆心之间的距离的几何关系求出 180 个方向角度，旋转中心坐标 （50.2775,56.1050）,相邻探测器间距离为 0.2775mm。 对于问题二和问题三， 本文根据中心切片定理， 建立起投影物体的密度与其 二维傅里叶变换后的函数的关系模型。在每一个确定角度下，分别使用直接反投 影法和滤波反投影法 ( , )|( cos ,sin )|( )f x yf rrg r ( )( , ) ( , )| gff x y f 将变换后的函数作傅立叶逆

3、变换就可以得到原物体的密度函数， 并进行插值 填充由 matlab 生成吸收率矩阵。 问题二的介质几何中心为 （43.1641,47.4609） ， 形状大致是椭圆形的嵌套（见图 5-9）。问题二的 10 个点吸收率（附件 4 坐标 点自上到下）分别为：0、0、2.0053、1.9714、1.9613、1.9644、0、0、0、0。 图片见本文 4.2.2 模型求解部分。问题三的 10 个点吸收率（附件 4 坐标点自上 到下）分别为：0、0、0.6479、0、0、4.7242、0、10.5331、0、0。图片轮廓 形状类似于组织细胞图，图片见本文 4.2.3 模型求解部分。 对于问题四，问题

4、1 中的相邻探测器间距离参数，是影响成像精度的重要参 数，距离越小精确度越高。旋转中心在托盘上的位置和探测器的个数，是影响成 像稳定性的重要参数，增加单位长度的探测器个数，并且调整旋转中心和探测器 中心相对位置，使得探测器每转一个角度，接收到信息探测器尽可能在中部（避 免下次转过一个角度时，介质的投影信息溢出光屏）。探测器绕旋转中心扫过的 方向角度尽量稍微大于 180，使信息采集的更加完备，从而增加了 ct 工作的 稳定性。 关键词：关键词： 余弦轨迹 滤波反投影法 中心切片定理,傅里叶变换 傅里叶级数分 析 2 一、问题重述一、问题重述 1.11.1 ct(computed tomograp

5、hy)可以在不破坏样品的情况下，利用样品对射线能 量的吸收特性对生物组织和工程材料的样品进行断层成像， 由此获取样品内部的 结构信息。一种典型的二维 ct 系统如图 1 所示，平行入射的 x 射线垂直于探测 器平面，每个探测器单元看成一个接收点，且等距排列。x 射线的发射器和探测 器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转 180 次。 对每一个 x 射线方向，在具有 512 个等距单元的探测器上测量经位置固定不 动的二维待检测介质吸收衰减后的射线能量， 并经过增益等处理后得到 180 组接 收信息。 ct 系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的

6、ct 系统进行参数标定， 即借助于已知结构的样品 （称为模板） 标定 ct 系统的参数， 并据此对未知结构的样品进行成像。 1.21.2 问题的提出问题的提出 1.在正方形托盘上放置两个均匀固体介质组成的标定模板， 其中每一点的数 值反映了该点的吸收强度， 这里称为 “吸收率” 。 请根据这一模板及其接收信息， 确定 ct 系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该 ct 系统使用的 x 射线的 180 个方向。 2.利用(1)中得到的标定参数，确定附件 3 中的未知介质在正方形托盘中的 位置、几何形状和吸收率等信息。另外，请具体给出图 3 所给的 10 个位置处的 吸收率，

7、相应的数据文件见附件 4。 3. 附件 5 是利用上述 ct 系统得到的另一个未知介质的接收信息。利用(1) 中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图 3 所给的 10 个位置处的吸收率。 4. 分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立 对应的标定模型，以改进标定精度和稳定性，并说明理由。 二、问题分析二、问题分析 2.12.1 问题一分析问题一分析 以正方形托盘左下角为坐标原点建立平面直角坐标系。 讨论光线方向角和 探测器与 x 轴夹角的关系，根据投影方向角、旋转中心的坐标、旋转中心到最 近的探测器的标号之间的关系，建立光路方程。光路方程被模板

8、的几何图形所 截长度正比于所给的吸收率。经过旋转，圆心的运动轨迹随角度的变化是余弦 图并且其直流分量 k 即为旋转中心所对应的探测器的编号， 对其曲线拟合和傅 3 里叶级数分析求出 k 值。联立两次 k 探测器接收的光线的直线方程，交点就是 旋转中心，投影数据分析求探测器间距、180 个方向。 2.22.2 问题二分析问题二分析 问题二是一个已知角度下的投影，还原原函数密度的问题。根据中心切片 定理： 函数密度在某一方向上的投影的一维傅里叶变换函数是原函数密度的二 维傅里叶变换函数在同一角度且过原点直线上的值。 据此定理结合附件数据求 解变换后的傅里叶函数在每一个角度下的频域函数，再取傅立叶逆

9、变换。 2.32.3 问题三分析问题三分析 问题三与问题二分析和求解类似， 并且都是建立在第一问的已知 180 次测量 角度上的， 由于上述处理方法非特殊取值的方法， 因此该方法具有普遍适用型， 此处不再赘述。 2.42.4 问题四的分析问题四的分析 模板对稳定性的影响在于模板在正方形托盘位置的分布， 由几何概型可知， 模板分布越集中于托盘中部，被探测器检测到的概率越大，稳定性越高。因此 本文设计出的模板模型是圆心位于正方形托盘重心处的圆。 三、模型假设三、模型假设 1.忽略空气对 x 射线的吸收率。 2.不考虑光的折射和反射。 3.假设所测量的数据均为正常测量的情况 ，不考虑测量过程的抖动因

10、素。 4.假设 ct 的滤波器是理想的。 四、符号说明四、符号说明 :光线的方向角，光线矢量与 y 轴正方向夹角,逆时针旋转为增大的正方向 :探测器和 x 轴正方向夹角 i:探测器编 （当探测器旋转到平行于正方形上方时， 从左向右以 0512 编号） k:距离旋转中心最近的探测器编号 q:吸收率 d:椭圆中心和圆心之间对应多个探测器的距离 l:相邻探测器间距离 4 y1,y2:包络线位置函数 max1 y、 max2 y：椭圆中心和圆心对应探测器的编号 n：测量次数（1n180） w：傅里叶展开函数的基频率 （a,b）：旋转中心坐标 1 f、 2 f：圆的 8mm 直径对应 28 个间距、29

11、 个间距的探测器，出现频率 ( , )f x y：未知物体的密度函数 ( )gr ：密度函数在方向上的投影函数 ( )g：密度函数在方向上的投影函数一维傅里叶变换函数 ( , )f ：( , )f x y的二维傅里叶变换函数 ( , ) w fx y：重建后的密度函数 ( )h w：滤波器的频域表达式 ( , ) wh f ：( , )f 滤波后的结果 ( )n：每一次测量角度的弧度制 12 tt、:存储相邻两次旋转的稳定因子矩阵 t l：程序中的矩阵的长度 n：值为 512，代表 512 个探测器 五五、模型的建立与求解、模型的建立与求解 5.15.1.1.1 问题问题 1 1 的模型建立的

12、模型建立 建立平面直角坐标系建立平面直角坐标系求解求解 和和经过经过 k k 探测器的光路探测器的光路直线方程直线方程（分三类讨论）（分三类讨论） （1）090 o 5 图 5-1（090 o旋转过程图） 由几何关系可知： = cos d l arccos d l 经过 k 探测器的光路直线方程：tan(tan ) kk yxbag （2）90180 oo 图 5-2（90180 oo旋转过程图） 同理可得： 90 o 6 d sin = l =arcsin90 d l o 经过 k 探测器的光路直线方程： tan(90 )tan(90 ) kk yxba oo g （3）180270 oo

13、图 5-3（180270 oo旋转过程图） 同理可得 180 o d sin = l =arcsin180 d l o max1max2 |dyylg 经过 k 探测器的光路直线方程：tan(180 )tan(180 ) kk yxba oo g 求求 k k 由几何关系可知，对于单一的几何图形来说，在任一确定角度上，通过椭圆中 心的弦和圆的圆心的弦是最长的，即对应的探测器接收数据最大。将附件 2 的数 7 据，对每一列数据使用牛顿迭代法找到极大值 max1 y、 max2 y极大值对应的探测器 编号即为穿过圆心的光线对应的探测器。横轴代表 180 次测量的转次（暂未知角 度），纵轴代表一次测

14、量时，每一个的传感器测量状态（状态：传感器有无吸收 率，1i512）。 图 5-4 圆与椭圆双旋转图 曲线拟合方法： 由于两条包络是围绕同一旋转中心得到的，所以它们的相位相等，将两条包 络线进行取平均得到 y。 图 5-5 椭圆和圆连线中点在接收器上的投影位置图 8 假使探测器不动，模板相对探测器旋转，圆心绕旋转中心的轨迹，分时投影到 探测器上的曲线是余弦函数曲线形状， 旋转后的结果图如下：（旋转动图见附件） 。 5-6 几何中心点围绕旋转中心旋转产生的余弦图像 对它利用傅里叶级数展开可以得到直流分量和余弦分量， 直流分量的大小就是 k。 101 cos wsin w n yaanbn 进行一

15、级傅立叶级数展开： 1101 coswsinwyaab 11 20 cos2wsin2wyaab 11 30 cos3wsin3wyaab 101 cos wsin w n yaanbn 其中 y1、y2、y3已知，a0、a1、b1是未知。 用上面求出的 值对应的每一次位置，计算出每一次经过 k 探测器的光路直线方程， 求两两直线的焦点，求出旋转中心。由平面几何分析可得，每一次经过 k 探测器的光路 直线方程为： tan(tan ) kk yxbag 联立任意两个角度下的经过 k 探测器的光路直线方程，计算出多组两直线的交点， 利用统计求取均值，所得就是旋转中心点。 由于圆在探测器上投影的永远

16、是直径，为定值，因此在每一个角度确定的情况下，通 过圆直径所对应的探测器的间距个数，求出l 12 88 2928 lff 9 5.1.25.1.2 问题问题 1 1 模型的求解模型的求解 根据数据按列迭代求取极大值， 从而确定过椭圆中心和圆心每一个 下对应的 d 的大小，最后利用 4.1.2 中三个情况下的 数学表达式，求出每一个角度，得 到下表（见附件“180 个”.xls）： n n12345678910 角度角度29.498431.574431.574432.892734.165734.786835.398437.180537.758438.8923 n n111213141516171

17、81920 角度角度39.999140.542841.612943.177243.688644.697346.177652.217852.661253.5405 n n21222324252627282930 角度角度53.976354.840955.269856.120956.963657.382058.625559.445459.852761.0649 n n31323334353637383940 角度角度61.865262.263163.448264.621165.396866.167667.315468.075368.831169.9577 n n414243444546474849

18、50 角度角度70.704371.447472.556273.658274.389475.118076.206477.289678.009279.0851 n n51525354555657585960 角度角度80.157081.225581.936082.999684.060789.957989.957989.957989.957989.9579 n n61626364656667686970 角度角度89.957989.957989.957989.957989.957989.957995.586196.292897.354698.0640 n n71727374757677787980

19、角度角度99.4865100.1999101.2731101.9908103.0709103.7936104.8820105.6106106.7084107.4438 n n81828384858687888990 角度角度108.5526109.6686110.4169111.5463112.3041113.0661114.2172114.9905115.7686116.9455 n n919293949596979899100 角度角度117.3405118.5342119.3375120.1473120.9637121.7872122.6180123.0364123.8791124.73

20、02 n n101102103104105106107108109110 角度角度125.5902126.4595126.8979127.3387128.2284134.3117135.8048136.8228137.3390138.3871 n n111112113114115116117118119120 角度角度139.4572140.0010141.1077142.2416143.4053144.6016145.2132146.4654147.1073148.4256 n n121122123124125126127128129130 角度角度149.1036150.5016151.2

21、237151.9627152.7201154.2965155.1195156.8472156.8472157.7584 n n131132133134135136137138139140 角度角度158.7065159.6967161.8311161.8311162.9947164.2414164.2414167.0810167.0810167.0810 n n141142143144145146147148149150 角度角度167.0810170.7370170.7370170.7370173.2637173.2637173.2637177.7557173.2637186.7363 n

22、n151152153154155156157158159160 角度角度182.2443186.7363182.2443186.7363182.2443186.7363189.2630186.7363189.2630189.2630 n n161162163164165166167168169170 角度角度191.2388191.2388192.9190192.9190194.4075195.7586197.0053197.0053199.2646199.2646 n n171172173174175176177178179180 角度角度199.2646201.2935202.241620

23、3.1528204.0313204.8805206.5025206.5025208.0373208.0382 对 y 傅里叶展开求解： 101 cos wsin w n yaanbn 求以上解方程组再进行近似可得： a0 =256.8 95%置信区间(254.7, 259) a1 = 87.92 95%置信区间(86.64, 89.19) b1 = -76.42 95%置信区间(-80.75, -72.08) w = 0.01718 95%置信区间(0.0169, 0.01746) 因为探测器为整数个，故取 a0=257。 10 其中 是直流分量。 所以 k=257，第 257 个探测器到旋转

24、中心距离最短 871814 0.2775mm 29852885 l 5.25.2.1.1 问题问题 2 2 的模型建立的模型建立 （1）直接反投影法： 由第一问可知每一次测量的角度。 因此问题二的求解是建立在第一问的结果上 的。根据中心切片定理，对于每一个方向（对应于附件表格的每一列）都有： ( , )|( cos ,sin )|( )f x yf rrg r ( )( )gg r f ( )( , ) ( , )| gff x y f （即中心切片定理） 学术实验研究表明： 1 ( , )( , )* w fx yf x y r 。左式傅里叶变换： 1 ( , )( , ) w ff 其中

25、1 r 为造成还原图像模糊的原因，为消除模糊的现象，需要在频域上串联一 个滤波器，滤波器的函数表达式为( )h w。 下面介绍第二种方法滤波反投影法 为输出的( , ) wh f 串联一个滤波器( )h w ( , )( , )( ) wh ffh w g 得： ( , )( , ) wh ff 使得达到更强的精准特性。 5.2.25.2.2 问题问题 2 2 模型的求解模型的求解 将附件 3 按列导入，应用 4.2.14.2.1 问题问题 2 2 模型建立模型建立的数学表达式，的数学表达式， 使用 matlab 软件运行后，生成一个 256256 吸收率矩阵，通过计算图形过 x 轴与 y 轴

26、投影的中点的直线的交点得到介质的几何中心坐标为 （43.1641,47.4609），得到这 10 个点的吸收率结果为： 11 横坐标纵坐标吸收率 10.000018.00000.0000 34.500025.00000.0000 43.500033.00002.0053 45.000075.50001.9714 48.500055.50001.9613 50.000075.50001.9644 56.000076.50000.0000 65.500037.00000.0000 79.500018.00000.0000 得到的图像为： 图 5-7 直接反投影结果图 图 5-8 滤波反投影结果图

27、12 附件 3 只有 180 个数据，由于角度少，拟合效果不好，我们采取了增加几 个合理角度来生成更加精准图， 并且新的图比图 4-8 更加美观， 新图如下图 4-9： 图 5-9 拟合之后的图 5 5.3.3 问题问题 3 3 的模型建立的模型建立与求解与求解 问题三的模型建立与求解过程与问题二的完全一致。 得到的拟合生成图为图 5-10 图 5-10 问题 3 的生成图 13 所得的吸收率分别为： 横坐标纵坐标吸收率 10.000018.00000.0000 34.500025.00000.0000 43.500033.00000.6479 45.000075.50000.0000 48.

28、500055.50004.7242 50.000075.50000.0000 56.000076.50000.0000 65.500037.000010.5331 79.500018.00000.0000 5 5.4.4 问题问题 4 4 模型的建立模型的建立 建立一个圆心和正方形重心重合的圆形模板，由于存在对称性，因此若 ct 是稳定的并且精准的，则每经过相同角度，探测器接收到的信息也具有对称性。 若在相邻探测器间距很小的情况下，探测器接收到的信息仍然均为对称性的，那 么由物理实验误差理论可知：该 ct 仪的绝对误差必定小于相邻探测器间距离的 大小。由此得出结论：缩小相邻探测器间的距离，是减

29、小误差提高精度的最有效 的方法。标定模板使用中心对称度较高的图形（圆）是最佳选择。 六六、模型的评价、模型的评价 6.1 模型的优点 （1）本文对问题 1 提供的模型简单易懂。 （2）问题 2 引入傅里叶和中心切片定理等内容，都属于现代技术，能很好的贴 合实际，应用性好。 6.2 模型的缺点 对于滤波法重塑原物体模板，因为存在衰减特性，现实不存在这样理想的滤 波器可供使用。 14 七、七、参考文献参考文献 【1】 李保磊, 张耀军. x 射线 ct 系统投影旋转中心的测量j. 光学精密工 程, 2011, 19(5):967-971. 【2】 廖小丽, 李永杰, 陈李容. ct 原理 课程双语

30、教学的探索与实践j. 实验科学与技术, 2012(6):266-268. 【3】 王岩岩, 吕振伟. 一种改进的滤波反投影图像重建方法j. 科协论坛, 2007(8):1-2. 【4】 黄敏. ct 成像和胶片误差的分析j. 江苏大学学报(医学版), 1996(zj):68-68. 15 附录 matlab 程序代码 问题问题 1 1： % y1=peaks1(:,1) ; % x2=peaks2(:,2) ; % y2=peaks2(:,1) ;%求两条曲线函数散点 % subplot(2,1,1),plot(x1,y1,r:); % hold on ; % subplot(2,1,2),p

31、lot(x2,y2,b-p); % p1=polyfit(x1,y1,4); % p2=polyfit(x2,y2,4);%绘图 % p1=p1 ; createfigure1(peaksy, peaksx); %粗略椭圆和圆在 接收器上投影图 ，生成过程参考 createfigure1。m 文件 % 傅里叶椭圆和圆连线中点在接收器上的投影模拟图生成过程 n=347; for i=1:n if peaks(i,2)=peaks(i+1,2) 16 i=i+2; n=n-1; end end end %舍弃不合理的局部极大值 n=347; for i=1:n if peaks(i,2)=peak

32、s(i+1,2) road(i)=(peaks(i+1,1)+peaks(i,1)/2; end if peaks(i,2)=peaks(i+1,2) end end %取椭圆和圆连线中点在接收器上的投影位置 j=1; for i=1:346 if road(i)=0 y1(j)=road(i); j=j+1; end end plot(y1); x=1:179; 17 fitresult, gof = createfit(x, y1); %生成傅里叶拟合椭圆和圆连线中点 在探测器上的投影位置图 cftool; % 出现对话框，点击 date，选 x,y1,选择 fourier 模拟 %以傅里

33、叶模型拟合，求出偏移值 a0 % general model fourier1: % f(x) = a0 + a1*cos(x*w) + b1*sin(x*w) %coefficients (with 95% confidence bounds): % a0 = 256.8 (254.7, 259) %a1 = 87.92 (86.64, 89.19) % b1 = -76.42 (-80.75, -72.08) % w = 0.01718 (0.0169, 0.01746) function fitresult, gof = createfit(x, y1) %createfit(x,y)

34、% create a fit. % % data for untitled fit 1 fit: % x input : x % y output: y % output: 18 % fitresult : a fit object representing the fit. % gof : structure with goodness-of fit info. % % 另请参阅 fit, cfit, sfit. % 由 matlab 于 16-sep-2017 11:36:43 自动生成 % fit: untitled fit 1. xdata, ydata = preparecurved

35、ata( x, y1); % set up fittype and options. ft = fittype( fourier1 ); opts = fitoptions( method, nonlinearleastsquares ); opts.display = off; opts.startpoint = 0 0 0 0.0352987938605595; % fit model to data. fitresult, gof = fit( xdata, ydata, ft, opts ); % plot fit with data. figure( name, untitled f

36、it 1 ); h = plot( fitresult, xdata, ydata ); legend( h, y vs. x, untitled fit 1, location, northeast ); % % 创建 axes % axes1 = axes(parent,figure); % hold(axes1,on); 19 % 创建 xlabel xlabel(转次); % 创建 ylabel ylabel(椭圆和圆连线中点在探测器位置); % box(axes1,on); % grid(axes1,on); % % 创建 legend % legend(axes1,show); %

37、 % 创建 arrow % annotation(figure1,arrow,0.903571428571429 0.958928571428571,. % 0.110904761904762 0.10952380952381); % % % 创建 arrow % annotation(figure1,arrow,0.13125 0.13125,. % 0.926195426195427 0.986486486486487); grid on function createfigure(y1) %createfigure(y1) % y1: y 数据的矢量 % 由 matlab 于 16-se

38、p-2017 15:55:40 自动生成 20 % 创建 figure figure1 = figure; % 创建 axes axes1 = axes(parent,figure1); hold(axes1,on); % 创建 plot plot(y1,displayname,y1); % 创建 xlabel xlabel(转次); % 创建 ylabel ylabel(椭圆和圆连线中点在探测器器上的投影位置); box(axes1,on); % 设置其余坐标轴属性 set(axes1,outerposition,0 0 1 1); % 创建 legend legend(axes1,show

39、); 问题问题 2 2： %phantom %randon clear all; %p = phantom(modified shepp-logan,159); %r=radon(p); %figure;imshow(r,); %figure; 21 %imshow(p,);title(); r=as2; % %lt = pow2(nextpow2(size(r,1)-1);% lt=360; p1 = zeros(lt,lt);% for i=1 : size(r,2) temp = imrotate( repmat(r(:,i),1,size(r,1),i-1,bilinear ); te

40、mp = temp(floor(size(temp,1)/2-lt/2)+1:floor(size(temp,1)/2+lt /2),floor(size(temp,2)/2-lt/2)+1:floor(size(temp,2)/2+lt/ 2); p1=p1+temp; end p1=p1/size(r,2); p1=rot90(p1); figure;imshow(p1,);title(); % n=180; % h=size(r,1); h=zeros(h*2-1),1); for i=0:h-1 if i=0 h(h-i)=1/4; elseif rem(i,2)=0 h(h-i)=0

41、; h(h+i)=0; else h(h-i)=-1/(i*pi)2; h(h+i)=-1/(i*pi)2; end end x=zeros(h,n); for i=1:n s=r(:,i); xx=conv(s,h); x(:,i)=xx(h:2*h-1); end 22 % p3=zeros(lt,lt); for i=1:lt for j=1:lt for k=1:180 theta= (k+30)/180*pi; t=(j-lt/2-0.5)*cos(theta)+(lt/2+0.5-i)*sin(theta)+(h+1)/2 ; t1=floor(t); t2=floor(t+1);

42、 p3(i,j)=p3(i,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k); end end end p3=pi/n*p3; figure;imshow(p3,);title(); 问题问题 3 3 代码代码： clear all; a = xlsread(a 题附件.xls,附件 1,a1:iv256); as1 = xlsread(a 题附件.xls,附件 2,a1:fx512); as2 = xlsread(a 题附件.xls,附件 3,a1:fx512); as3 = xlsread(a 题附件.xls,附件 4,a1:b10); as4 =xlsread(a 题附件.xls,附件 5,a1:fx512); r=as4;% 读