搜索
CT系统参数标定及成像研究摘要本文通过建立数学模型，给出了通过分析大量数据标定CT参数的方法，得到了通过数据借助模型和相关的算法及函数生成图像的方法；以及其逆向运算方法。针对问题一，主要是对附件一、附件二中数据的分析和处理，标定CT系统的各个参数。根据附件一中模板的吸收率，对样品做连续化处理得到相应模板的方程，建立坐标系，关键词：连续化、参数标定、最小二乘原理、radon算法、成像一、问题重述CT(ComputedTomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。一种典型的二维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个位置示意图二、问题分析1、问题一该问题主要为数据分析处理类问题。首先对附件一进行分析，通过对模板每一点吸收率连续化处理分别得到椭圆和圆的方程。分析附件二中数据，并用MATLAB画出对应矩阵的图像，每一列不为0的数据个数为𝑀𝑖，反映第i次扫描有介质区域的长度，对比模板，反映探测器与模板的相对位置，𝑀𝑖最大时探测器平行于y轴，𝑀𝑖最时探测器平行于x轴；模板由均匀固体介质组成，则经过增益后的数据与实际数据成比例，附件中每个单元格中数据即为增益后的截模板弦长的量值。选取表格中数据量较大处的数据与相应的实际数据用近似最小二乘定理算得误差最小时的比例系数μ三、模型假设1、假设CT系统的旋转中心位于探测器中垂线上2、假设探测器上每一点发出光强相同，且均匀分布在整个探测器上3、忽略光的衍射四、符号系统𝑀𝑖：第i列不为0的数据个数(i=0,1,2,…,180)x轴：与椭圆短半轴重合的数轴y轴：过椭圆中心且与长半轴重合的数轴𝑥0：旋转中心横坐标𝑦0：旋转中心纵坐标d：探测器单元之间的距离𝐵𝑖𝑗：第i次旋转第j个探测器单元射线的截距𝑃𝑖𝑗：表格中第i列第j行的数据(i=1,2…，180，j=1,2,…,512)𝐿𝑖𝑗：探测器第i次扫描时第j个单元探测器扫描到的介质弦长和(i=1,2…，180，j=1,2,…,512)μ：问题一中弦长与衰减后射线能量的比例系数𝑘𝑖：射线斜率(i=1,2,…,180)θi：旋转第i次与X正半轴的夹角(i=0,1,2,…,180)𝑛𝑗：第n个探测器单元（j=1,2,3，…，512）五、模型建立1、问题一该问题主要为数据分析处理类问题。首先对附件一进行分析，通过对模板每一点吸收率连续化处理分别得到椭圆和圆的方程与理论一致分别为：𝑥2152+𝑦2402=1，(𝑥−45)2+𝑦2=16分析附件二中数据，并用MATLAB画出对应矩阵的图像（如图一），根据图像，探测器逆时针旋转，则探测器处于平行于y轴时，沿x负半轴方向分别为第1个到第512个探测单元。𝑀𝑖反映第i次扫描有介质区域的长度，反映探测器与模板的相对位置，对比模板实际位置，𝑀𝑖最大时探测器平行于y轴，𝑀𝑖最时探测器平行于x轴；模板由均匀固体介质组成，则经过增益后的数据与实际数据成比例，附件中每个单元格中数据即为增益后的截模板弦长的量值，即𝑃𝑖𝑗=μ𝐿𝑖𝑗。1）计算μ首先选取𝑀𝑖最大时的六列数据，可以有效减小系统误差，取每一列数据数值最大的即探测器位于平行于y轴的位置，最大值表示椭圆短轴和圆直径吸收衰减后的射线能量经增益处理的量值，取六个方向平均值L1=67.3419，对应的实际值为38；同理选取𝑀𝑖最小时的三组数据，此时探测器位于平行于x轴的位置，两段不为0数据中的最大值分别表示椭圆长半轴和圆直径吸收衰减后的射线能量增益后的量值，取三个方向平均值分别得到L2=141.6893667,L3=14.17846667，对应的实际数值为80和8;为得到更加精确的结果，对这三个数据用近似最小二乘法算误差最小时取得的比例系数，最终算得μ=1.7713，LINGO程序见附录一。2）计算d选取数据时首先考虑数据量大的位置，可以有效避免系统误差，提高精度。计算探测器单元间距d时，选取探测器𝑀𝑖最大时的六组𝑀𝑖，即探测器位于平行于y轴的位置，对取平均值𝑀𝑖得到椭圆长轴所对应的探测器单元个数，并对应长轴实际长度算出探测器单元间距：𝑑1=80289=0.276816609选取Mi最小时的三组数据，即探测器平行于x轴的位置，由于这几组数据中间部分有探测器单元没有扫描到介质区域，所以不能直接选取选取Mi计算，分别选取模板中椭圆右侧端点与圆左侧端点、椭圆左侧端点与圆右侧端点及圆和椭圆两短轴端点间对应的探测器单元个数，分别对应实际的数据得到三组探测器单元间距：𝑑2=2694=0.276595744𝑑3=64∗3694=0.27665706𝑑4=38∗3412=0.276699029对于得到的四组探测器单元间距，为了减小系统误差，求平均值保留四位小数得到d=0.2767，详细数据见附件表一。3）计算旋转中心坐标计算探测器旋转中心坐标时，利用已经计算出的d值和两条不相交的直线确定一个点得到旋转中心坐标。取探测器平行于y轴和平行于x轴两个位置的数据，这些数据中量值最大的对应的即为椭圆的中点，在这两个位置将椭圆中心即坐标系原点与旋转中心之间的探测器单元数目差值分别确定，找到模板和探测器系统的相对位置，代入d值，分别求得纵坐标和横坐标。探测器平行于y轴时椭圆中心对应的探测器单元：𝑛𝑦=(243+239+236+233+229+226)6⁄=233.3333𝑦0=（256-𝑛𝑦）*d=6.2719探测器平行于x轴时椭圆中心对应的探测器单元：𝑛𝑥=(222+223+223)3⁄=22.6667𝑥0=（𝑛𝑥-256）*d=-9.2233详细数据见附件表二，旋转中心相对坐标原点的位置示意图如图二。图二旋转中心坐标示意图4）CT系统射线的180个方向为确定探测器旋转的180个方向，设过旋转中心的探测器单元的射线方程为𝑦−𝑦0=𝑘(𝑥−𝑥0)，则任意一个探测器单元在任意旋转次数时的射线方程一般式为：𝑦−𝑦0=𝑘(𝑥−𝑥0)+(256−𝑛𝑗)∗𝑑𝑐𝑜𝑠θi(n=1,2…,512,i=1,2…,180)根据表格中第一列和最后一列数据，即探测器初始和最终位置，第一列中间部分数据为0的即为与椭圆和圆相切的两条射线间的探测器单元数∆𝑛，列出两条切线方程分别为：{𝑦=𝑡𝑎𝑛θ1∗𝑥+𝐵1𝑗𝑦=𝑡𝑎𝑛θ1∗𝑥+𝐵1𝑗+∆𝑛∗𝑑𝑐𝑜𝑠θ1初始位置与椭圆和圆相切射线方程分别与椭圆和圆方程联立{𝑦=𝑡𝑎𝑛θ1∗𝑥+𝐵1𝑗𝑥2152+𝑦2402=1{𝑦=𝑡𝑎𝑛θ1∗𝑥+𝐵1𝑗+∆𝑛∗𝑑𝑐𝑜𝑠θ1(𝑥−45)2+𝑦2=16相切则判别式∆=0可得到225∗𝑡𝑎𝑛θ12−𝐵1𝑗2+1600=02009∗𝑡𝑎𝑛θ12−90∗𝑡𝑎𝑛θ1∗(𝐵1𝑗+∆𝑛∗𝑑𝑐𝑜𝑠θ1)+(𝐵1𝑗+∆𝑛∗𝑑𝑐𝑜𝑠θ1)2−16=0利用探测器单元间距d,用MATLAB求得初始射线的角度，程序见附录一，0并将弧度制转化为角度制，得到初始角度θ1=119.7305，同理可得θ180=298.7164。相减得总旋转角度θ=θ180−θ1=178.9859°.将射线一般式方程分别与椭圆方程和圆的方程联立求解{𝑦−𝑦0=𝑘(𝑥−𝑥0)+(256−𝑛𝑗)∗𝑑𝑐𝑜𝑠𝜃𝑖𝑥2152+𝑦2402=1{𝑦−𝑦0=𝑘(𝑥−𝑥0)+(256−𝑛𝑗)∗𝑑𝑐𝑜𝑠𝜃𝑖(𝑥−45)2+𝑦2=16利用韦达定理和弦长公式算出含有k和𝜃𝑖的弦长表达式，其中弦长公式为L=√1+𝑘2√(𝑥1+𝑥2)2−4∗𝑥1∗𝑥2，k=√1−𝑐𝑜𝑠𝜃𝑖2𝑐𝑜𝑠𝜃𝑖.与椭圆相交弦长：𝐿𝑖𝑗1=√1+𝑘2(9∗𝑘2+64)⁄*√(𝑦0+(256−𝑛𝑗)∗𝑑𝑐𝑜𝑠𝜃𝑖)2−18𝑘2𝑥0−4∗(9∗𝑘2+64)(9∗(𝑦0+(256−𝑛𝑗)∗𝑑𝑐𝑜𝑠𝜃𝑖)2−1600∗9)与圆相交的弦长：Lij2=2∗√16−(𝑘∗(45−𝑥0)+(256−𝑛𝑗)∗𝑑𝑐𝑜𝑠𝜃𝑖)21+𝑘2利用已经计算出的μ值，将附件二中数据处理后代入到这512*180个方程中，每一个i值对应的512个方程写成一个方程组，由最小二乘定理使得该方程组误差最小时取得相应的180个θ值，用mathematica求解，一共可解得512组解，分别取平均值得到最后的结果。方程组为𝐿𝑖𝑗=𝐿𝑖𝑗1+𝐿𝑖𝑗2(𝑖=1,2,…,180;𝑗=1,2,…,512)由于数据量巨大，计算机不能够正确运算出结果，计算时选择数据量最大的长半轴处进行求解最为精确，又因为椭圆弦长和圆弦长不可分离，实际计算导致部分方程组不能求出可行解故构造分段函数求解，以是否通过小圆区域作为临界条件，将180组氛围三段计算，临界时必有射线与小圆相切，从而解得临界射线斜率分别为-0.0418和-0.1908，程序见附录二，三组方程分别为{𝐿𝑖𝑗=𝐿𝑖𝑗1𝐿𝑖𝑗=𝐿𝑖𝑗1+𝐿𝑖𝑗2𝐿𝑖𝑗=𝐿𝑖𝑗1180次旋转分别与x轴正方向夹角见表（1）119.7146121.0580121.6095122.6918123.7178124.6815125.6760126.6708127.6657128.6609129.6562130.6517131.6475132.6431133.6388134.7845135.6306136.6265137.6224138.6183139.6142140.6101141.6059142.6017143.5973144.5929145.5883146.5837147.5788148.5738149.5687150.4628151.5576152.5516153.5454154.5389155.5