行粗调；
[0156] 对影像产生扭曲需说明的是，在安装CCD相机时，除了步骤D中提到的底座不平问 题，而且会在水平或竖直方向产生一定的倾斜角度，这便会使采集到的离子影像变成椭圆。 本发明中，处理影像变椭问题采取的是影像压缩的方式（将椭圆影像的长轴进行压缩）。另 外，由于本步骤中仅用一个参数对影像的扭曲进行调节，故调节结果会显得略微粗糙，因此 被称为粗调。
[0157] 影像变椭有两种：影像的横向比纵向较宽或影像的纵向比横向较宽。若影像的横 向较宽，则应将其横向进行压缩，即采取横向压缩的方式；若影像的纵向较宽，则应将其纵 向进行压缩，即采取纵向压缩的方式。
[0158] 由于纵向压缩的原理与思想与横向压缩完全一样，故以下的陈述内容仅以横向压 缩为例，纵向压缩便不再赘述。
[0159] 图7为图1所示处理方法中影像压缩粗调步骤的实施例流程图。请参照图7,以横 向压缩为例，该粗调步骤E包括：
[0160] 步骤701 :设定横向压缩系数coeff ;
[0161] 该横向压缩系数coeff的取值范围介于0~1之间。coeff取值越接近0,影像压 缩程度越小；coeff取值越接近1，影像压缩程度越大。
[0162]以下各步骤中，运用数学变换，将旋转0^角度后的三维离子速度聚焦影像所对应 的矩阵转换为一个新矩阵，该新矩阵所对应的影像便为压缩处理后的影像，其比原影像稍 圆。
[0163]步骤702:将i、j都赋初值为0;
[0164] 步骤703:设置循环步，i = i+1作为外循环；
[0165] 步骤704:将j = j+1作为内循环；
[0166] 步骤705:表达式i-(i_xQ)*coeff表示i在压缩后的新位置，由于i-(i-xQ)*coeff 可能为小数，故需要对其进行取整，此处采用小数点后第一位四舍五入的取整方法。先引入 变量 k，再将 i-(i_xQ)*coeff 取整后的值即 round(i-(i_xQ)*coeff)赋给 k;
[0167] 步骤706:将（i，j)处的矩阵数值I(i，j)赋给"（k，j)，其中"为mXn的0 矩阵，实现了（i，j)在横向的压缩；
[0168] 步骤707 :判断是否j < n，若是，执行步骤704 ;否则，内循环结束，执行步骤708 ;
[0169] 步骤708:判断是否i < m，若是，执行步骤703;否则，外循环结束，执行步骤709;
[0170] 步骤709:画圆，将圆作为影像压缩过程中的对比参照物。观察压缩后影像与所画 圆形之间在所有方向的差异是否都在五个像素点之内，若是，则说明压缩后影像与所画圆 形之间契合程度较好，压缩参数输入合适，粗调完成；若否，则说明压缩后影像与所画圆形 之间契合程度不好，则需更换压缩参数重复以上压缩步骤702~708,直至压缩后影像与所 画圆形之间在所有方向的差异都被控制在五个像素点之内为止。
[0171] 步骤G:对完成粗调后的三维离子速度聚焦影像进行扭曲的微调，使影像变得更 圆；
[0172] 由于飞行时间质谱中的电场分布不会是理想状态，会产生一定的扭曲，那么，采集 到的三维离子速度聚焦影像不仅会因为CCD的摆放角度产生扭曲，而且会因为电场分布产 生其他更复杂的扭曲。对于更复杂的扭曲，则需要更精细的调节。
[0173] 本发明中，对三维离子速度聚焦影像进行微调，需通过以下一系列的数学公式来 实现：
[0174] i = ij+di (s, t) (1)
[0175] di = a1t+a2t2+a3t3 (2)
[0176] j = ji+dj (s, t) (3)
[0177] dj = b1s2+b2s3+b3ts 2 (4)
[0178] sQ) = 2(i「x0)/m (5)
[0179] t(j!) = 2(j1-y〇)/n (6)
[0180] 其中，X(l，y(l为矩阵中影像中心对应的行和列，m，n为整个矩阵的行数和列数，i，j 为原影像所对应的矩阵的某行和某列，ii，1为处理后影像所对应的矩阵的某行和某列，s， t为公式变换的中间变量。系数ap a2, a3, bp b2, b3为可调参数，调节这6个系数便可对整 个矩阵进行定量变换。与粗调相比，该步骤中用六个参数对影像的扭曲进行调节，其调节过 程更加精细，调节结果更加完美，因此被称为微调。
[0181] 在实现整个变换的过程中，其中一个较大的难点便是由原影像对应矩阵的（i，j) 计算出微调后对应的（ip 。因为matlab计算能力有限，不能对较复杂的公式实施变换 和推导，故需要先进行人为推导，将公式变换成多项式形式，再利用matlab程序平台上现 有的求多项式根的函数，便可计算出（ii，i)，其具体推导过程如下：
[0182] 将⑵式代入⑴式得：
[0183] i = i1+a1t+a2t2+a3t 3 (7)
[0184]由（5)式得：
[0185] it =(^) ^ + ^0 (8)
[0186] 将⑶式代入（7)式得：
[0187]
[0188] 再将⑷式代入（3)式得：
[0189] j = j1+b1s2+b2s3+b 3ts2 (10)
[0190]由（6)式得：
[0191]
[0192]
[0193]
[0194] 此时，联合（9) (12)两式便可得到一个二元三次方程组，理论上是可以直接解出 变量s和t的值，但是由于matlab计算能力有限，无法直接解出方程组的解。因此需要人 为变换方程组的形式，推导出一个matlab可以直接获得变量s和t解的形式。经过多番考 虑，在matlab程序中有一个求解多项式根的rootsO函数可直接被调用。那么，为了可以 利用roots ()函数来进行求解，则需要将（9) (12)两式组合成的方程组进行变量代换，得出 一个单变量多项式。
[0195] 由（9)式得：
[0199] 以上（14)式便是一个关于t且最高项为9次方的多项式方程，那么利用roots () 函数便能直接求出t的解。当求出t的值之后，将t值代入（6)式便可求出的值。
[0200] 之后将t值代入（9)式，可求出相应的s值。再将s值代入（5)式则可求出对应 的L值。
[0201] 至此，整个推导过程结束，（ip i)值被计算出来。
[0202] 图8为图1所示处理方法中影像定量变换微调步骤的实施例流程图。请参照图8， 具体而言，该步骤G包括：
[0203] 步骤801 :输入参数&1，a2, a3, bp b2, b3的设定值；
[0204] 步骤802 :定义一个与I (i，j)维数相同的0矩阵A，即A = 1*0 ;
[0205] 步骤803 :将i、j都赋初值为0 ;
[0206] 步骤804 :设置循环步，i = i+1作为外循环；
[0207] 步骤805 :将j = j+1作为内循环；
[0208] 步骤806 :利用下述方法计算出（i，j)在微调后对应的（ip i):
[0209] (1)计算么的具体过程如下：
[0210] 首先，对以下方程组进行变量代换，将其转换成关于t的多项式方程，再由t的多 项式方程求解出t的值：
[0211]
[0212] m，n分别是整个矩阵的行数和列数，s、t是中间变量，（xQ，yQ)为影像圆环中心的 坐标，a" a2, a3, bp b2, b3为输入的调节参数。
[0213] 而后，利用t值由下式计算j1:
[0214] t(j!) = 2(j1-y〇)/n
[0215] (2)计算il的具体过程如下：
[0216] 首先，利用t值由下式计算出s的值：
[0217]
[0218] 再利用求出的s值由下式计算i1:
[0219] sQ) = 2(i「x0)/m
[0220] 其中，上述计算公式的推导过程在前已经详细说明，此处不再重述。
[0221] 步骤807 :由于ip i可能不为整数，故需对i i进行取整，且取整后的值分别赋 给i'，j'，本步骤中采取向下取整的方式。
[0222] 步骤808 :在进行取权重计算时会用到ip i分别与i'，j'的差值，故需将差值 分别赋为S i'，6 j';
[0223] 步骤809 :再次使用取权重的方法将影像所对应的矩阵实施变换。
[0224] 为了使微调变换后的影像各点更加连贯，此处再次使用取权重的方式进行操作， 即
[0225] (I)A(ir )=A(i'，j，)+I(i, j)*(l-8 ir )*(l-8 )；
[0226] (II)A(i, +l)=A(ir +1)+I(i, j)*(l-8 ir )*8j,；
[0227] (IlDAa，+1，j，）= A(i，+1，j，）+I(i，j)* S i，*(1-S j，）；
[0228] (iyMU' +1，j，+1) = A(i' +1，j，+1)+I(i，j)* S i，* S j，。
[0229] 步骤810 :判断是否j < n，若是，执行步骤805 ;否则，内循环结束，执行步骤811 ;
[0230] 步骤811 :判断是否i < m，若是，执行步骤804 ;否则，外循环结束，则产生的矩阵 便是微调后的影像所对应的矩阵。
[0231] 步骤H:对微调后的三维离子速度聚焦影像，计算其速度分布，依据该速度分布的 分辨率判断是否要继续进行微调，如果是，则需调整输入的 ai，a2, a3, bp b2,匕的设定值，重 新执行步骤G ;否则，执行步骤I ;
[0232] 前述步骤C为计算三维离子速度聚焦影像的角分布，计算角分布是：将矩阵I (r， 9)某固定角度处位于所求半径范围内的强度相加并乘以积分因子。那么对于速度分布的 计算，便是将矩阵I(r，0)某固定半径处位于所求角度范围内的强度相加并乘以积分因子 即可得出某固定半径处的强度分布。
[0233] 图9为图1所示处理方法中计算影像速度分布步骤的详细流程图。请参照图9,具 体分析，该步骤中部分处理方法和过程与步骤C是一样的。
[0234]如图所示，其步骤901，902，903，904，905，907，908与图4所示一样，只有步骤906， 909,910,911有所区别。由于计算速度分布与角分布存在着积分因子的差异，因此角度 ~9max范围内的速度分布计算如下：将矩阵1(:，9 中的所有矩阵元乘以积 分因子r2|sin0 |(步骤906)，再将矩阵的列进行累加，便得到角度0min~0 _范围内速 度分布"（r)(步骤909)。
[0235] 具体而言，该步骤H具体包括：
[0236] 步骤901 :运用坐标转换公式：x = xQ+rcos 0，y = yQ+rsin 0将极坐标r，0与直 角坐标x，y联系起来；
[0237] 步骤902 ;将r，0赋初值为0 ;
[0238] 步骤903 :设置循环步，r = r+1作为外循环；
[0239] 步骤904 :将0 = 0+1作为内循环；
[0240] 步骤905 :运用多点取权重的方法将直角坐标（i，j)中的矩阵I (i，j)转换至极坐 标（r，0 )中的矩阵I(r，0 )，即将（r，0 )周围四个点的强度的权重加起来赋给I(r，0 )， 其中：
[0241] (I) I (i (r，0)，j (r，0 ))的权重为（1-8 j (r，0 ))*(1_ 8 i (r，0 ));
[0242] (Il)I(i(r，0)，j(r，0 )+l)的权重为 S j(r，0 )*(1-S i(r，0 ));
[0243] (Ill)I(i(r，0)+l，j(r，0))的权重为（1-S j(r，0))*Si(r，0);
[0244] (IV) I (i (r，0 )+l，j (r，0 )+l)的权重为 8 j (r，0 )* 8 i (r，0 );
[0245] 故而，I(r，0) = I(i(r，0)，j(r，0))*[(l-8j(r，0))*(l-8i(r， Q))]+I(i(r，0)，j(r，0)+l)*[8j(r，0)*(l-8i(r，0))]+I(i(r，0)+l，j(r， 0))*[(l-6j(r，0))*8i(r，0)]+I(i(r，0)+l，j(r，9)+l)*[8j(r，0)*8i(r，0)]
[0246] 步骤 906 :令 I ' （r，9 ) = I (r，9 ) *r21 sin 9 | ;
[0247] 步骤907 :判断是否0 < 360,如果是，执行步骤904 ;否则，内循环结束，执行步骤 9