CT系统参数标定数学模型

CT系统参数标定数学模型摘要本文通过对体模系统和投影点坐标的分析，提出了较为充分合理的假设，并对数据进行了拟合处理。先后建立了两个模型（理想状态下的模型，只考虑载物台和探测器等一些因素误差的模型）来求解CT系统参数标定的问题。在对问题一的分析中，只考虑在理想情况下：即载物台水平和探测器无偏转，并不考虑机械误差的情况下各个参数之间的关系，我们以探测器左下角为（0，0），右方向为x轴正方向，上方向为y轴正方向，垂直于xoy平面的走直线为z轴建立如图4所示的空间直角坐标系。运用几何知识作图，通过建立相似图形的等比例方程来研究确定几个系统参数之间的关系（X射线源的位置参数、载物中心线的位置参数和探测器的位置参数）。最后讨论了引入载物台中心轴倾斜、载物转台不均匀转动和数据不确定度对误差的影响。在分析问题二时，考虑到所给椭圆点集可能有偏差，于是首先运用MATLAB软件对数据进行处理拟合，得到更为确信的数据并由此求出修正后的“椭圆”方程。从拟合出来的椭圆方程，我们可以得知椭圆的中心坐标还有其大小，并因此看出我们拟合出来的椭圆并不和我们的理想模型相符合。所以我们想到了有机械误差存在。由于同时考虑多个误差因素将使模型计算变得极为复杂，对此，我们假设探测器无偏转，只有载物台倾斜，建立模型二，并对探测器上的两个椭圆进行旋转修正。易知修正前后两个椭圆的大小并没有改变，只是相对位置发生了变化，发射源A在探测器上的投影坐标也没有改变。对此，我们以A点为旋转中心对其进行旋转修正并根据拟合后的数据求出两个类椭圆中心连线的斜率。从而求出载物台倾斜角ａ=arctan0.00210563最后在问题三，对模型三的机械误差进行了分析计算。机械误差主要来源为探测器上数据像素误差、载物台的转轴偏斜，圆心距离的误差，并对其分析了原因。关键词CT系统载物台转轴偏斜数据修正拟合机械误差参数标定问题重述CT自发明以来，被公认为自伦琴发现X射线后是放射领域的最重要的发明之一。工业CT一般采用扫描物件旋转和面探测器的结构，工作示意图如图1示：图SEQ图\*ARABIC1工业CT系统工作原理示意图扫描物件围绕某一固定转轴旋转，每隔一定角度采集一张图像，然后根据采集的图像采用3D图像重建算法即可将原始3D物件重建出来。在实际中，由X射线源发出X射线，经过扫描物体衰减后照射在探测器上，探测器根据接收到的光子数的计数实现光电转换，从而形成灰度图像。一般平板探测器大小为3000×2000像素，每个像素为0.127mm。在CT系统安装过程中往往存在机械误差，而这些误差对于物件重建的准确性往往是致关重要的，实际中就需要对安装好的CT系统进行参数标定，系统参数主要包括X射线源的位置、载物中心线的位置和探测器的位置参数等。传统的方法是采用尺子测量法，对于CT系统参数标定精度要求是不够的（通常CT重建过程中系统参数的机械误差不允许超过一个像素）。目前，一般通过实验的方法实现参数标定。所谓实验的方法标定参数就是通过扫描已知参数的体模，分析投影数据，估计系统参数的方法。针对工业CT要求解决如下问题：（1）建立合适的坐标系，正确描述CT系统的各种参数和机械误差，并建模分析这些参数的关系和可能的机械误差。（2）通常采用轴承钢球作为CT参数标定体模。因为它具有各向投影一致，边缘清晰等优点，在里采用两个钢球实现CT的参数标定实验（如图3所示）。两个钢球置于有机玻璃管中（X射线容易透过有机玻璃，容易后期球心获取的图像处理），钢球直径为8mm±0.0008mm，两球心距离约为100mm±1mm。将有机玻璃管固定在旋转载物台上，载物台携带有机玻璃管以均匀速度旋转，与此同时X射线源发出射线，探测器采集数据（旋转一周，等间距采集180张）。附件1给出了根据180张采集图像提取出的球心投影坐标。两个钢球球心投影的轨迹是两个椭圆（如图4所示）。试根据实验数据估计该CT系统的参数值，即给出CT系统的标定。（3）在参数标定过程中可能存在多种可能的机械误差，试就你的参数标定可能的误差进行分析。图SEQ图\*ARABIC2双球体模型示意图图SEQ图\*ARABIC3双球投影轨迹（双球轨迹坐标见附录1）问题分析题目要求我们求出工业CT系统的各项参数，这些参数包括以下六个：=1\*GB3①放射源到转盘的转轴的水平距离L1；=2\*GB3②转盘的转轴到探测器平面的水平距离L2；=3\*GB3③放射源到探测器平面的水平距离L3；=4\*GB3④放射源到xoz平面的垂直距离Ya；（平面xoz具体位置见问题（1）解答中的坐标系）=5\*GB3⑤放射源与yoz平面的水平距离Xa；=6\*GB3⑥转轴到yoz平面的水平距离d。如图所示：为了使参数更容易确定，我们选取探测器作为参照物，这样各个参数就有了参考的标准。问题一要求我们求出系统各个参数之间的关系，还有根据我们所列的关系式进行误差分析。本题误差主要来源于工业CT系统安装过程。除了装置本身的由于在装系统时，可能会导致物台中心轴倾斜、载物转台不均匀转动，探测器前后俯仰等，都会导致机械误差的存在。我们在计算时会考虑这些因素,并建立相应的模型。所以，在标定时应该给出其机械误差大小。可以根据体模的投影数据，利用几何知识来求解主要包括X射线源的位置、载物中心线的位置和探测器的位置并求出他们的关系。问题二则以题目给出的360个坐标为基础，运用matlab软件来拟合两个椭圆，考虑到这些点的坐标本身就有一点的误差，经过一定的修改和筛选后，进行多次拟合进而求出椭圆的方程。求出椭圆的方程后，再利用在问题（1）中各个参数的关系，就可以求出各个参数的值。根据分析，我们知道了拟合出来的两个椭圆的中心连线并不垂直于X轴，我们可以利用绕点旋转的方法对其进行修正。问题三中的误差分析，本题中考虑的机械误差来源有两圆心距离误差，载物台的中心转轴倾斜，这些根据问题二中的数据，我们已经对它进行了修正。模型的假设与符号说明模型假设假设光子波动不明显，即忽略光的衍射现象；将两个钢球的球心看作一个质点；假设实验所收集的数据能客观反映实际情况；玻璃管在随着转盘转动的过程中，始终垂直固定于转盘；载物台在转动过程中，玻璃管与两个钢球始终保持相对静止。符号说明及名词解释符号符号含义数据类型L放射源到转轴的水平距离待求参数L转轴到探测器平面的水平距离待求参数L放射源到探测器平面的水平距离待求参数X放射源到yoz平面的距离待求参数Y放射源到xoz平面的垂直距离待求参数d转轴到yoz平面的距离待求参数r有机玻璃管和小球的轨迹圆半径中间值a上椭圆的短轴中间值a下椭圆的短轴中间值b两椭圆长轴的平均值中间值m两椭圆最高点间的距离中间值k两椭圆最低点间的距离中间值t放射源在探测器上的投影与上椭圆最低点的距离中间值h两钢球球心间的距离已知量模型建立及问题求解问题一建立合适的坐标系，正确描述CT系统的各种参数和机械误差，并建模分析这些参数的关系和可能的机械误差。将探测器所在的平面定为参照物，结合题目已给出的二维坐标系xoy，在探测器平面左下方原点（0，0）引一条Z轴，形成三维坐标系，如下图示：依照题意，我们首先得求出载物台的转盘半径r，为了求出转盘的半径，我们先作了这样的一个实验，我们以探测器左下角为（0，0），右方向为x轴正方向，上方向为y轴正方向，垂直于xoy平面的走直线为z轴建立如图1所示的空间直角坐标系。，再让发射源与玻璃管中的下边小球处于同一个小平面上进行实验。在实验过程中，我们能在探测器上得到玻璃管中上面小球所形成的一条直线为2b，而下面小球的投影则形成一个类椭圆，我们可以通过探测器上的数据拟合得到椭圆的方程，从而得到我们想要椭圆的长轴长和短轴长还有中心坐标，进而可以计算出有机玻璃管和小球的轨迹圆半径r。图4为了简单分析，我们分别画出图1的正视图和俯视图如下：俯视图正视图根据正视图和俯视图，我们利用相似三角形和勾股定理得到了以下关系式：L1+L1-rrb=联立（1）（2）（3）可解得：r=b*h2*m*h1根据以上的式子，我们可以根据实验得出的数据求出r的值，然而对于相同的这个装置来说，以下的实验中有机玻璃管和小球的轨迹圆半径r是不会改变的，还有因为其俯视图中的2b的值是发射源垂直于xoy平面发射出来的投影的长度，而在其他情况下是求不出那个值，为了计算方便，在以下的运算中，我们将有机玻璃管和小球的轨迹圆半径r作为一个常数a进行计算，即令r=a（a则为由上面（4）求出的数值）接着，我们按照题目的要求保持发射源A与载物台和探测器的水平距离不变，即相对位置不变，而在竖直方向上移动发射源A，使其与两个小球都不在同一个水平面进行以下实验，其三维图如下：为了方便分析，我们作出它的正视图如下：根据其正视图，利用相似三角形定理和勾股定理可以得到以下关系式：L1+rL1-ryB-r=a（其中a为一个常数）（8）根据CT系统的工作原理，即令扫描物件围绕某一固定转轴旋转，并每隔一定角度采集一张图像，然后根据采集的图像采用3D图像重建算法将原始3D物件重建出来。再根据重建出来的图像来标定我们的系统参数，由此可知，凡是在重建图像所需要的数据误差都会影响我们的参数标定。因此，根据上方我们所列出来的关系式可以看出此过程中存在的误差有以下这几种：载物台没有正常匀速转动，引起重建图像不连续扫描物件的放置发生倾斜或在运动过程产生抖动，引起重建图像不规则载物台转动角度发生偏离或载物台的转轴发生倾斜，引起重建图像的数据不好处理实验中两个钢球的半径大小有一定的误差，使得有机玻璃管和小球的轨迹圆半径r存在误差从而影响参数标定。问题二根据照题目给出的数据，我们用MATLAB软件进行拟合，绘制出上下两个椭圆的图形（算法详见附录2），由图形反映拟合情况，并解出椭圆一般方程的系数p1、p2数据图形拟合上方椭圆的中心坐标为（1421.7199，1812.5604，0）长半轴长为1351.2899，短半轴长为172.8553（单位为像素）下方椭圆的中心坐标为（1427.2216，188.4994，0）长半轴长为1342.6651，短半轴长为上方椭圆的方程为：(x-1421.7199)下方椭圆的方程为：(x-1427.2216)由此，我们可以得出B,D,C,E各点的纵坐标分别为y将其坐标代入（5）（6）（7）（8）式可得：LLLya=811.83225，x即我们可以得出发射源A的坐标为（1423.832，811.83225，11.58469ａ）综上所述，我们可以对给定CT系统参数进行标定（单位为像素）：X射线源的位置与载物台中心线的距离：LX射线源的位置与探测器的距离：L若以坐标和方程的形式给出，则表示为：发射源A的坐标为（1423.832，811.83225，11.58469ａ）载物台中心线的方程x=探测器的位置：x-o-y平面模型的进一步分析与修正在上述的分析过程中，我们发现经过拟合出来的两个椭圆的中心连线并不垂直于X轴，与实际中我们预测分析的有所不同。由于影响系统的机械误差有很多，在此，我们只认为这是因为载物台在转动过程中所绕的转轴倾斜产生偏角所导致的，与其它因素无关。所以，我们对其作了修正如下（左为修正前，右为修正后）：由问题二我们已经求出了发射源A的坐标为（1423.832，811.83225，11.58469），由于此过程中只在xoy平面上修正，所以修正后发射源的坐标是不会改变的。对此，我们以A点在xoy平面上的投影A’为中心，线段OBOC绕着中心A’旋转到线段OBOC垂直于X轴，此时，上下两个椭圆的中心坐标分别为OB（1423.832，1812.56265，0上方椭圆的方程：(x-下方椭圆的方程：(x-此时我们还可以求出修正的角度ａ，其修正的图形在xoy平面的角度ａ如下图所示：以上方椭圆的中心坐标为标准，则修正前的角度ａ有如下关系：tan代入数据可以求出：tanａ即修正的角度的大小为：ａ=arctan经过了这样的修正，我们可以更加准确的利用探测器上的数据进行拟合图像，使得重建的图像更加具有准确性，我们能够更好的进行标定CT系统的参数。问题三除了问题二中所提到的载物台的转轴偏斜引起误差之外，还存在以下两种误差可能引起参数标定不准确。首先，如果载物台没有正常匀速转动，那么扫描物件的旋转速度会受到影响，探测器所采集到的具体数据可能会发生断点或重点，如下图的断点情况示意图：其次，扫描物件的放置发生倾斜或在运动过程产生抖动会导致探测器所采集到的图像数据发生偏差，影响了Matlab的拟合效果以及相关数据的准确性。如图12所示，很明显地，投影得到的图形发生了扭曲，形成机械误差，影响到了参数的标定。模型评价与改进优点：对体模的投影点坐标进行数据的修正，使得计算的结果更加准确；考虑到的误差比较合理，有较大的实际意义运用MATLAB软件进行数据运算和图表拟合，使得计算更加简单准确。缺点：机械误差存在的因素考虑还不够全面，且仅限于理论上的分析，有待实践过程中的进一步检验与改进在数据误差处理时仅用理论知识进行处理，实际操作有些难度。参考文献[1]百度文库/[2]中国知网/kns50/index.aspx[3]谢照鸿，范正森等人编著，《数学建模技术》，中国水利水电出版社，2003年9月[4]周永正，詹棠森，方成鸿，邱望仁，《数学建模》，同济大学出版社，2011年1月附录附录1双球模型投影在探测器上的坐标上方球心坐标下方球心坐标x轴y轴x轴y轴2762179227372142754178627272172745178127152202735177527012242723176926872272710176426712312695175826532342679175326352372661174726152402643174225942432623173725722462602173225492492580172725252522557172225002552532171724742572507171324472602481170824192622454170423902652425170023602672396169623302692367169222992712336168822672742305168422352752273168122022772240167721682792207167421342812173167120992832138166820642842103166520282852068166319922872032166019552881995165819182891959165518812901921165318432911884165118052921846165017672931807164817282941769164616902941730164516512951691164416122951652164315732961613164215332961573164114942961534164014552961494164014152961454164013762961415163913372961375163912972961336164012582951296164012192951256164011802941217164111412931178164111032931139164210642921100164310262911061164498829010231646951289985164791428894716498772869101651841285873165280528383616547692828001657734280764165970027872916616662776941664632275660166760027362616705682715931673536268561167650626652916794762644991683447261468168641825943916903912564101694364254383169833925135617023142483301706291245305171126824228117152472392591720227236237172520723321617301902301971735173226179174015822316217451442191461750131216132175612021211917611102091071767102205971773952028917789019882178486194761790841907217968318770180285183691808881797018149217573182199171771827107167841833117164921839129160102184514315611318521581521271858176148143186419514516018702171411791877240137201188326513422418892911302491895320127276190035112330419063831203351912417117367191745311440119234901104381928529108475193357010551419386121025551942656100598194770197642195174795687195579593734195984491782196389389832196694487882197099686933197210498598619751103841039197711578310931979121182114819811267811203198313228112591984137881131519851434811371198514898114271986154581148419861600821540198516558315961985171083165219841764841707198318178617611981187087181519801922891868197819739119211975202392197219732071952022197021199720721967216599212019632210102216619602253104221219562295107225519522336110229719482374113233919432412116237819382447119241519342480123245019292512126248419232542130251619182570133254719132596137257519072621140260219012643144262618962664148264918902683152267018842699156268918782714159270618722727163272218662739167273518592748171274618532756175275618472761179276418412765183277018352768186277418282768190277618222767194277718162765198277618102760201277318042754205276917982747209注：探测器左下角为（0，0），右方向为x轴正方向，上方向为y轴正方向附录2用matlab软件求椭圆方程的程序和输出结果：求上椭圆方程（先将上椭圆的180个点的坐标导进，记作x1，y1，用矩阵表示）F=@(p,x)p(1)*x(:,1).^2+p(2)*x(:,2).^2+p(3)*x(:,1)+p(4)*x(:,2)+p(5);p0=[0.0050.0050.0050.0050.005];x=[x1,y1]warningoffp=nlinfit(x,zeros(size(x,1),1),F,p0)plot(x(:,1),x(:,2),'ro');holdon;xmin=min(x(:,1));xmax=max(x(:,1));ymin=min(x(:,2));ymax=max(x(:,2));ezplot(@(x,y)F(p,[x,y]),[-1+xmin,1+xmax,-1+ymin,1+ymax]);formatlongpA=sqrt(p(3)^2/(4*p(1)^2)+p(4)^2/(4*p(1)*p(2))-p(5)/p(1))B=sqrt(p(3)^2/(4*p(1)*p(2))+p(4)^2/(4*p(2)^2)-p(5)/p(2))x0=p(3)/(2*p(1))y0=p(4)/(2*p(2))C=A^2D=B^2输出结果：p=1.0e-022*0.00000.0000-0.0000-0.00010.1249p=1.0e-022*0.000000000621500.00000003798134-0.00000176718735-0.000137686937760.12490432520128A=1.351289908892398e+003B=1.728