CN105321206A - 一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法 - Google Patents

Abstract

本发明公开了一种中子层析成像旋转轴线偏摆角的误差补偿方法，本方法是基于平行束中子射线的投影几何对称性原理，将样品360度旋转范围内的投影数据进行叠加获得具有唯一对称轴的合成投影图像，合成投影图像的对称轴即为样品旋转轴线的投影线；进一步统计出合成投影图像的像素值梯度角直方图，该直方图分布为一偶函数。求取该偶函数的对称中心，该对称中心位置对应的梯度角即为旋转轴线的偏摆角。将测量所得的偏摆角作为修正参数引入到重建算法中，即可实现对旋转轴线偏摆角的误差补偿，从而保证样品各重建断层清晰度的一致，有效提高重建精度。

Description

一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法

技术领域

本发明涉及一种中子层析成像旋转轴线偏摆角的误差补偿方法，可用于中子照相领域中的无损检测及工业计算机断层扫描。

背景技术

中子照相技术在检测含氢材料、重金属组件结构、原子序数相近的不同元素、放射性材料等方面有着X射线无法具备的优势，因而在军事、航空、航天、考古等领域有着重要的应用。作为中子照相的一个重要分支——中子CT(ComputedTomography，计算机断层成像术)也日益得到重视和快速的发展，并在检测样品内部结构及质量状态方面发挥了重要的作用。

图1为传统中子层析成像系统的层析扫描原理图，中子射束1为一组平行线，被检测样品3在旋转平台2的带动下绕旋转轴线4旋转360度，探测器5采集被检测样品3在不同旋转位置的投影图像，然后利用成像系统中的图像重建模块反演出样品的断层图像。由图1可知，图像重建模块所采用的图像重建方法是基于样品旋转轴投影线41要严格平行于探测器坐标系Y_d轴。在进行图像重建时，各个断层的反投影点仅是探测器坐标系X_d轴上的变量，且该反投影点的计算是以样品旋转轴投影线41为基准。由此可知，当旋转轴投影线41平行于探测器坐标系Y_d轴时，该投影线在X_d轴上的坐标为一常数，因此各个断层的反投影地址的基准也均等于一个常数，此种情形下进行断层图像重建时，采用统一的标准重建算法，即可确保重建出的各断层图像清晰度保持一致。

但是，在中子层析成像系统的实际搭建中，由于机械部件的制造和安装误差，导致样品旋转轴投影线41不平行于探测器坐标系Y_d轴，即样品旋转轴投影线41相对于探测器坐标系Y_d轴产生了一定的偏角，如图2所示。此时样品旋转轴投影线41在X_d轴上的坐标就不等于常数，从而导致各个断层的反投影点的基准也不为常数。在图像重建时如果仍按照常数进行计算，就会导致Y_d轴方向上所重建断层的清晰度分布不均，影响三维重构模型的精度和特征信息测量的可信度。

该问题解决的关键就是精确标定出样品旋转轴投影线41与探测器坐标系X_d轴的夹角，本发明中将该角度称为旋转轴线偏摆角α。为此，本发明设计出了一种针对中子层析成像系统样品旋转轴线偏摆角的自动测量方法。将测量所得的偏摆角作为修正参数引入到重建算法中，即可实现对旋转轴线偏摆角的误差补偿，从而保证样品各重建断层清晰度的一致，达到提高三维重构分辨率的目的。

发明内容

本发明方法是基于平行束中子射线的投影几何对称性原理，将样品360度旋转范围内的投影数据Pⁿ进行阈值分割以消除背景信息的干扰，然后将分割后的所有投影图像进行叠加获得具有唯一对称轴的合成投影图像P^ave。根据平行束投影的几何属性可得合成投影图像P^ave的对称轴即为样品旋转轴投影线，该投影线与探测器坐标系X_d轴的夹角即为旋转轴线偏摆角α；进一步统计出合成投影图像P^ave的像素梯度角直方图，该直方图分布为一偶函数F(β)；求取函数F(β)的对称中心位置对应的梯度角β₀，β₀即为旋转轴线的偏摆角α。将测量所得的偏摆角α作为修正参数引入到重建算法中，即可实现对旋转轴线偏摆角的误差补偿，从而保证样品各重建断层清晰度的一致，达到提高三维重构分辨率的目的。

本发明是一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法，其特征在于包括有下列步骤：

步骤一、样品投影信息的获取；

将投影图Pⁿ分割后的投影图记为分割投影图分割投影图与投影图Pⁿ的宽度和高度是相等的，第n个分割投影图内任意一点的投影数据记为那么Ostu阈值分割可表示为：

$p_{Th}^n(x_i,y_j)=\left\{\begin{array}{cc}{p}^n(x_i,y_j)& p^n(x_i,y_j)\&GreaterEqual;TValue\\ 0& p^n(x_i,y_j)<TValue\end{array}\right.$

表示第n个分割投影图 $P_{Th}^n$ 的第i行第j列的采样点投影值；

pⁿ(x_i,y_j)表示探测器采集到的第n个投影图Pⁿ的第i行第j列的采样点投影值；

TValue为基于类间方差最小准则获取的最优分割阈值；

所述 $P^n={\left[\begin{array}{cccccc}{p}^n(x_1,y_1)& p^n(x_1,y_2)& ...& p^n(x_1,y_j)& ...& p^n(x_1,y_w)\\ p^n(x_2,y_1)& p^n(x_2,y_2)& ...& p^n(x_2,y_j)& ...& p^n(x_2,y_w)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ p^n(x_i,y_1)& p^n(x_i,y_2)& ...& p^n(x_i,y_j)& ...& p^n(x_i,y_w)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ p^n(x_h,y_1)& p^n(x_h,y_2)& ...& p^n(x_h,y_j)& ...& p^n(x_h,y_w)\end{array}\right]}_{h\&times;w};$

步骤二、合成投影图的获取；

将所有分割投影图 $P_{Th}^n$ 叠加，则获得另一幅二维投影图像，该二维投影图像即为上述的合成投影图像P^ave，P^ave内的任意一个投影值p(x_i,y_j)的计算式为：

$p(x_i,y_j)=\frac{1}{M}\underset{n=1}{\overset{M}{\&Sigma;}}{p}_{Th}^n(x_i,y_j)$

其中n＝1,2,3......M，M表示样品在360度旋转范围内的投影总数目。

所述 $P^{ave}={\left[\begin{array}{cccccc}p(x_1,y_1)& p(x_1,y_2)& ...& p(x_1,y_j)& ...& p(x_1,y_w)\\ p(x_2,y_1)& p(x_2,y_2)& ...& p(x_2,y_j)& ...& p(x_2,y_w)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ p(x_i,y_1)& p(x_i,y_2)& ...& p(x_i,y_j)& ...& p(x_i,y_w)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ p(x_h,y_1)& p(x_h,y_2)& ...& p(x_j,y_j)& ...& p(x_h,y_w)\end{array}\right]}_{h\&times;w};$

步骤三、获取梯度角及梯度角直方图；

求取合成投影P^ave的所有像素点的梯度角，P^ave内任意投影值p(x_i,y_j)梯度角β(x_i,y_j)的计算公式如下：

$\mathrm{\&beta;}(x_i,y_j)=atan\frac{{\&dtri;}_{x}p(x_i,y_j)}{{\&dtri;}_{y}p(x_i,y_j)}$

$\begin{array}{c}{\&dtri;}_{x}p(x_i,y_j)=p\&lsqb;p(x_i+1,y_j+1)+(x_i,y_j+1)+p(x_i-1,y_j+1)\&rsqb;\\ -\&lsqb;p(x_i+1,y_j-1)+(x_i,y_j-1)+p(x_i-1,y_j-1)\&rsqb;\end{array}$

$\begin{array}{c}{\&dtri;}_{y}p(x_i,y_j)=\&lsqb;p(x_i-1,y_j-1)+p(x_i-1,y_j)+p(x_i-1,y_j+1)\&rsqb;\\ -\&lsqb;p(x_i+1,y_j-1)+p(x_i+1,y_j)+p(x_i+1,y_j+1)\&rsqb;\end{array}$

其中，表示任意投影值p(x_i,y_j)在X_d轴方向的梯度值，表示任意投影值p(x_i,y_j)在Y_d轴方向的梯度值。β表示任意像素的梯度角，x_i为探测器坐标系X_d轴上的变量，y_j为探测器坐标系Y_d轴上的变量。

统计所有像素点的梯度角的直方图，梯度角直方图的计算式如下：

$F\left(\mathrm{\&beta;}\right)=\frac{n\left(\mathrm{\&beta;}\right)}{w\&times;h}$

其中F(β)表示直方图函数，该函数自变量为像素点的梯度角β，w表示合成投影图像P^ave的宽度，h表示合成投影图像P^ave的高度，n(β)表示梯度角等于β的像素点的个数。

步骤四、获取旋转轴线偏摆角；

求取梯度角直方图函数F(β)的中心点坐标β₀，该坐标即为旋转轴投影线41与X_d轴的夹角α，计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&beta;}}_0=\frac{1}{180}\underset{\mathrm{\&beta;}=0}{\overset{180}{\&Sigma;}}\mathrm{\&beta;}\&times;F\left(\mathrm{\&beta;}\right)$

步骤五、旋转轴线偏摆角误差补偿；

将计算出的旋转轴投影线41与X_d轴的夹角α作为修正参数引入到中子层析成像系统的重建算法中，即可实现对旋转轴线偏摆角的误差补偿。

本发明标定方法的优点在于：

①本发明解决了由于机械部件的制造和安装误差造成的样品转轴偏摆的问题，在保证样品各重建断层清晰度一致的同时，降低了重建图像像质对扫描平台制造和安装精度的要求，降低了生产和装配成本。

②本方法基于平行束中子射线的投影几何对称性原理，直接利用样品的原始投影信息即可反求出旋转轴线的偏摆角，不需要制作专门的校准模体和进行单独的扫描来实现误差补偿。

③本发明将旋转轴线偏摆角的测量转换为求取合成投影图像对称轴线倾角的问题。由于对称轴线的唯一性，从而保证了测量结果的唯一性，因此在实际实施中不易出现假解。

④本发明方法适用于平行束扫描模式下任何具有旋转扫描功能的CT系统，且实现中无需对CT扫描系统进行任何改动，即可嵌入现有CT重建软件中作为辅助升级模块，有效提高了CT系统的适应能力。

附图说明

图1是标准中子层析成像原理图。

图2是样品旋转轴线发生偏摆时的中子层析成像原理图。

图3是合成投影图像示意图。

图4是实际的合成投影图像。

图5是本发明合成投影图像的梯度角直方图。

图6是旋转轴线偏摆角误差未补偿时的重建断层图像。

图7是采用本发明方法后的旋转轴线偏摆角误差补偿后的重建断层图像。

1.中子射束	2.旋转平台	3.被检测样品
			4.旋转轴线	41.旋转轴投影线	5.平板探测器

具体实施方式

下面将结合附图和实施例对本发明做进一步的详细说明。

如图1所示，样品旋转轴线4虚拟不可见，其在探测器平面的投影41亦不可见。当旋转轴线投影41平行于探测器坐标系Y_d轴时，所述旋转轴线投影41与X_d轴的夹角为90度，且在探测器坐标系X_dO_dY_d中的方程为x＝x₀，x为探测器坐标系X_d轴上的变量，x₀为旋转轴投影线41与探测器坐标系X_d轴的交点坐标。在中子层析成像系统的实际搭建中，由于机械部件的制造和安装误差，导致样品旋转轴投影线41不平行于探测器坐标系Y_d轴，即样品旋转轴投影线41相对于探测器坐标系Y_d轴产生了一定的偏角。此时，旋转轴投影线41与X_d轴的夹角就不等于90度，如图2所示。本发明的发明点就是精确测量出旋转轴投影线41与X_d轴的夹角α，以实现对旋转轴线偏摆角的误差补偿。

被检测样品3在旋转平台2的带动下旋转360度，旋转平台2步进旋转，步进角步长一般取0.5°～1°。样品每转过一个角步长，探测器5采集到被检测样品3的一幅投影图，该投影图记为Pⁿ，即被检测样品的第n个投影，其中n＝1,2,3......M，M表示被检测样品在360度旋转范围内的投影总数目，一般取720～1200。可以看出，Pⁿ即为一个二维图像矩阵，可表示为：

$P^n={\left[\begin{array}{cccccc}{p}^n(x_1,y_1)& p^n(x_1,y_2)& ...& p^n(x_1,y_j)& ...& p^n(x_1,y_w)\\ p^n(x_2,y_1)& p^n(x_2,y_2)& ...& p^n(x_2,y_j)& ...& p^n(x_2,y_w)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ p^n(x_i,y_1)& p^n(x_i,y_2)& ...& p^n(x_i,y_j)& ...& p^n(x_i,y_w)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ p^n(x_h,y_1)& p^n(x_h,y_2)& ...& p^n(x_h,y_j)& ...& p^n(x_h,y_w)\end{array}\right]}_{h\&times;w}---\left(1\right)$

在式(1)中被检测样品第n个投影图Pⁿ内任意一点的投影数据记为pⁿ(x_i,y_j)，其中x_i为探测器坐标系X_d上的变量，i＝1,2,3...h，h表示投影图Pⁿ的高度；y_j为探测器坐标系Y_d轴上的变量，j＝1,2,3...w，w表示投影图Pⁿ的宽度。

pⁿ(x₁,y₁)表示探测器采集到的第n个投影图Pⁿ的第1行第1列的采样点投影值；

pⁿ(x₁,y₂)表示探测器采集到的第n个投影图Pⁿ的第1行第2列的采样点投影值；

pⁿ(x₁,y_j)表示探测器采集到的第n个投影图Pⁿ的第1行第j列的采样点投影值；

pⁿ(x₁,y_w)表示探测器采集到的第n个投影图Pⁿ的第1行第w列的采样点投影值；

pⁿ(x₂,y₁)表示探测器采集到的第n个投影图Pⁿ的第2行第1列的采样点投影值；

pⁿ(x₂,y₂)表示探测器采集到的第n个投影图Pⁿ的第2行第2列的采样点投影值；

pⁿ(x₂,y_j)表示探测器采集到的第n个投影图Pⁿ的第2行第j列的采样点投影值；

pⁿ(x₂,y_w)表示探测器采集到的第n个投影图Pⁿ的第2行第w列的采样点投影值；

pⁿ(x_i,y₁)表示探测器采集到的第n个投影图Pⁿ的第i行第1列的采样点投影值；

pⁿ(x_i,y₂)表示探测器采集到的第n个投影图Pⁿ的第i行第2列的采样点投影值；

pⁿ(x_i,y_j)表示探测器采集到的第n个投影图Pⁿ的第i行第j列的采样点投影值；

pⁿ(x_i,y_w)表示探测器采集到的第n个投影图Pⁿ的第i行第w列的采样点投影值；

pⁿ(x_h,y₁)表示探测器采集到的第n个投影图Pⁿ的第h行第1列的采样点投影值；

pⁿ(x_h,y₂)表示探测器采集到的第n个投影图Pⁿ的第h行第2列的采样点投影值；

pⁿ(x_h,y_j)表示探测器采集到的第n个投影图Pⁿ的第h行第j列的采样点投影值；

pⁿ(x_h,y_w)表示探测器采集到的第n个投影图Pⁿ的第h行第w列的采样点投影值。

在得到被检测样品360度旋转范围内的M幅投影图后，为了提取出每幅投影图Pⁿ中样品的完整投影信息，需要对投影图Pⁿ进行分割，以去除背景信息的干扰。本发明中所采用的分割方法为Ostu阈值分割法。所述的Ostu阈值分割法的原理是基于类间方差最小准则获取最优分割阈值，然后利用阈值将原图像分成前景、背景两个部分。将投影图Pⁿ分割后的投影图记为表示第n个分割投影图。可以看出，分割投影图与投影图Pⁿ的宽度和高度相等，第n个分割投影图内任意一点的投影数据记为那么Ostu阈值分割可表示为：

$p_{Th}^n(x_i,y_j)=\left\{\begin{array}{cc}{p}^n(x_i,y_j)& p^n(x_i,y_j)\&GreaterEqual;TValue\\ 0& p^n(x_i,y_j)<TValue\end{array}\right.---\left(2\right)$

表示第n个分割投影图的第i行第j列的采样点投影值；

pⁿ(x_i,y_j)表示探测器采集到的第n个投影图Pⁿ的第i行第j列的采样点投影值；

TValue为基于类间方差最小准则获取的最优分割阈值。

在本发明中，通过分析中子层析扫描的几何原理可知，将所有的分割投影图叠加，可获得另一幅二维投影图像，这里将该投影图像称为合成投影图像P^ave。可以看出，合成投影图像P^ave仍为一个二维图像矩阵，记为：

$P^{ave}={\left[\begin{array}{cccccc}p(x_1,y_1)& p(x_1,y_2)& ...& p(x_1,y_j)& ...& p(x_1,y_w)\\ p(x_2,y_1)& p(x_2,y_2)& ...& p(x_2,y_j)& ...& p(x_2,y_w)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ p(x_i,y_1)& p(x_i,y_2)& ...& p(x_i,y_j)& ...& p(x_i,y_w)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& & \&CenterDot;\\ p(x_h,y_1)& p(x_h,y_2)& ...& p(x_h,y_j)& ...& p(x_h,y_w)\end{array}\right]}_{h\&times;w}---\left(3\right)$

式(3)中合成投影图像P^ave内任意一点的投影数据记为p(x_i,y_j)，其中x_i为探测器坐标系X_d上的变量，i＝1,2,3...h，h表示投影图P^ave的高度；y_j为探测器坐标系Y_d轴上的变量，j＝1,2,3...w，w表示投影图P^ave的宽度。

p(x₁,y₁)表示合成投影图像P^ave的第1行第1列的采样点投影值；

p(x₁,y₂)表示合成投影图像P^ave的第1行第2列的采样点投影值；

p(x₁,y_j)表示合成投影图像P^ave的第1行第j列的采样点投影值；

p(x₁,y_w)表示合成投影图像P^ave的第1行第w列的采样点投影值；

p(x₂,y₁)表示合成投影图像P^ave的第2行第1列的采样点投影值；

p(x₂,y₂)表示合成投影图像P^ave的第2行第2列的采样点投影值；

p(x₂,y_j)表示合成投影图像P^ave的第2行第j列的采样点投影值；

p(x₂,y_w)表示合成投影图像P^ave的第2行第w列的采样点投影值；

p(x_i,y₁)表示合成投影图像P^ave的第i行第1列的采样点投影值；

p(x_i,y₂)表示合成投影图像P^ave的第i行第2列的采样点投影值；

p(x_i,y_j)表示合成投影图像P^ave的第i行第j列的采样点投影值；

p(x_i,y_w)表示合成投影图像P^ave的第i行第w列的采样点投影值；

p(x_h,y₁)表示合成投影图像P^ave的第h行第1列的采样点投影值；

p(x_h,y₂)表示合成投影图像P^ave的第h行第2列的采样点投影值；

p(x_h,y_j)表示合成投影图像P^ave的第h行第j列的采样点投影值；

p(x_h,y_w)表示合成投影图像P^ave的第h行第w列的采样点投影值。

在本发明中，合成投影图像P^ave具备如下属性：合成投影图像存在唯一的一个对称轴，该对称轴线即为旋转轴投影线41。由此可知，测量旋转轴投影线41与X_d轴的夹角α即就是寻找合成投影图像P^ave的对称轴。

由上述对称性可知，如图3所示，合成投影图像P^ave中的任意一个像素点记为像素点A，与所述像素点A对称的像素点记为投影像素点B，则存在唯一的投影像素点B与所述像素点A关于旋转轴投影线41对称。像素点A的投影值记为p(x_A,y_A)，梯度角记为β^A，x_A表示像素点A在探测器坐标系X_d轴上的坐标，y_A表示像素点A在探测器坐标系Y_d轴上的坐标；投影像素点B的投影值记为p(x_B,y_B)，梯度角记为β^B，x_B表示投影像素点B在探测器坐标系X_d轴上的坐标，y_B表示投影像素点B在探测器坐标系Y_d轴上的坐标。则存在如下关系：

$p(x_A,y_A)=p(x_B,y_B),\mathrm{\&alpha;}=\frac{{\mathrm{\&beta;}}^A+{\mathrm{\&beta;}}^B}{2}---\left(4\right)$

进一步计算出合成投影图像P^ave中所有像素点的梯度角β，可以看出，β^A和β^B是β的两个取值。统计出梯度角的直方图，那么该直方图为梯度角β的函数，记为F(β)，且F(β)的分布为一偶函数，其对称点的位置所对应的梯度角记为β₀。则有，β₀就是旋转轴投影线41与X_d轴的夹角α。

获得旋转轴投影线41与X_d轴的夹角α后，即可将该角度作为修正参数引入到中子层析成像系统的重建算法中，实现对旋转轴线偏摆角的误差补偿，从而保证各重建断层清晰度的一致。

本发明是一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法，该方法包括有下列实施步骤：

步骤一、样品投影信息的获取；

获取样品在360度旋转范围内的投影图，即得到M幅投影图Pⁿ。为了提取出每幅投影图Pⁿ中样品的完整投影信息，采用Ostu阈值分割法对每幅投影图Pⁿ进行分割，以去除背景信息的干扰。

本发明中所述的Ostu阈值分割法是一种基于类间方差最小准则获取最优分割阈值，然后利用阈值将原图像分成前景、背景两个部分。将投影图Pⁿ分割后的投影图记为分割投影图分割投影图与投影图Pⁿ的宽度和高度是相等的，第n个分割投影图内任意一点的投影数据记为那么Ostu阈值分割可表示为：

$p_{Th}^n(x_i,y_j)=\left\{\begin{array}{cc}{p}^n(x_i,y_j)& p^n(x_i,y_j)\&GreaterEqual;TValue\\ 0& p^n(x_i,y_j)<TValue\end{array}\right.$

表示第n个分割投影图的第i行第j列的采样点投影值；

pⁿ(x_i,y_j)表示探测器采集到的第n个投影图Pⁿ的第i行第j列的采样点投影值；

TValue为基于类间方差最小准则获取的最优分割阈值。

步骤二、合成投影图的获取；

将所有分割投影图叠加，则获得另一幅二维投影图像，该二维投影图像即为上述的合成投影图像P^ave，P^ave内的任意一个投影值p(x_i,y_j)的计算式为：

$p(x_i,y_j)=\frac{1}{M}\underset{n=1}{\overset{M}{\&Sigma;}}{p}_{Th}^n(x_i,y_j)---\left(5\right)$

其中n＝1,2,3......M，M表示样品在360度旋转范围内的投影总数目，一般取720～1200。

步骤三、获取梯度角及梯度角直方图；

求取合成投影P^ave的所有像素点的梯度角，P^ave内任意投影值p(x_i,y_j)梯度角β(x_i,y_j)的计算公式如下：

$\mathrm{\&beta;}(x_i,y_j)=atan\frac{{\&dtri;}_{x}p(x_i,y_j)}{{\&dtri;}_{y}p(x_i,y_j)}---\left(6\right)$

$\begin{array}{c}{\&dtri;}_{x}p(x_i,y_j)=p\&lsqb;p(x_i+1,y_j+1)+(x_i,y_j+1)+p(x_i-1,y_j+1)\&rsqb;\\ -\&lsqb;p(x_i+1,y_j-1)+(x_i,y_j-1)+p(x_i-1,y_j-1)\&rsqb;\end{array}---\left(7\right)$

$\begin{array}{c}{\&dtri;}_{y}p(x_i,y_j)=\&lsqb;p(x_i-1,y_j-1)+p(x_i-1,y_j)+p(x_i-1,y_j+1)\&rsqb;\\ -\&lsqb;p(x_i+1,y_j-1)+p(x_i+1,y_j)+p(x_i+1,y_j+1)\&rsqb;\end{array}---\left(8\right)$

其中，表示任意投影值p(x_i,y_j)在X_d轴方向的梯度值，表示任意投影值p(x_i,y_j)在Y_d轴方向的梯度值。β表示任意像素的梯度角，x_i为探测器坐标系X_d轴上的变量，y_j为探测器坐标系Y_d轴上的变量。

统计所有像素点的梯度角的直方图，梯度角直方图的计算式如下：

$F\left(\mathrm{\&beta;}\right)=\frac{n\left(\mathrm{\&beta;}\right)}{w\&times;h}---\left(9\right)$

其中F(β)表示直方图函数，该函数自变量为像素点的梯度角β，w表示合成投影图像P^ave的宽度，h表示合成投影图像P^ave的高度，n(β)表示梯度角等于β的像素点的个数。

步骤四、获取旋转轴线偏摆角；

求取梯度角直方图函数F(β)的中心点坐标β₀，该坐标即为旋转轴投影线41与X_d轴的夹角α，计算公式如下：

$\mathrm{\&alpha;}={\mathrm{\&beta;}}_0=\frac{1}{180}\underset{\mathrm{\&beta;}=0}{\overset{180}{\&Sigma;}}\mathrm{\&beta;}\&times;F\left(\mathrm{\&beta;}\right)---\left(10\right)$

步骤五、旋转轴线偏摆角误差补偿；

将计算出的旋转轴投影线41与X_d轴的夹角α作为修正参数引入到中子层析成像系统的重建算法中，即可实现对旋转轴线偏摆角的误差补偿。

实施案例

实验采用的扫描装置参数如下：

(1)中子源：冷中子源

(2)面阵探测器：成像尺寸2048×2048，像元尺寸为0.2mm。

实验步骤：

(1)将一被检测样品放置在扫描平台上，启动扫描平台旋转以带动被检测样品步进旋转360°，步进角为0.5°，共采集720幅投影图。对所有投影图进行Ostu阈值分割，得到720幅分割投影图；然后将所有分割投影图叠加合成一个二维矩阵，该矩阵即为合成投影图像P^ave，如图4所示。

(2)根据本发明方法中的获取梯度角及梯度角直方图步骤，求取合成投影图像P^ave的所有像素点的梯度角，并统计出梯度角的直方图F(β)，如图5所示。可以看出，F(β)具有很好的对称性，则有F(β)对称中心位置对应的梯度角即为旋转轴线的偏摆角。

(3)根据获取旋转轴线偏摆角步骤求取梯度直方图函数的中心点坐标β₀，该坐标即为旋转轴投影线41与X_d轴的夹角α，该案例中α的实测值为87.7°。

(4)利用原始的投影进行重建，即对旋转轴线的偏摆角误差不进行补偿，默认其等于90度，重建结果如图6所示，可以看出，重建断层出现了明显的重影，降低了图像的空间分辨率。将旋转轴线偏摆角α＝87.7°代入重建软件的预处理校正模块，利用校正后的投影数据进行重建，结果如7所示。可以看出，图像的重影得到了很好的消除，图像清晰度得到了明显的提升，证明了本发明方法的有效性。

Claims (5)

1.一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法，其特征在于包括有下列步骤：

步骤一、样品投影信息的获取；

将投影图Pⁿ分割后的投影图记为分割投影图分割投影图与投影图Pⁿ的宽度和高度是相等的，第n个分割投影图内任意一点的投影数据记为那么Ostu阈值分割可表示为：

表示第n个分割投影图的第i行第j列的采样点投影值；

pⁿ(x_i,y_j)表示探测器采集到的第n个投影图Pⁿ的第i行第j列的采样点投影值；

TValue为基于类间方差最小准则获取的最优分割阈值；

所述

步骤二、合成投影图的获取；

将所有分割投影图叠加，则获得另一幅二维投影图像，该二维投影图像即为上述的合成投影图像P^ave，P^ave内的任意一个投影值p(x_i,y_j)的计算式为：

其中n＝1,2,3......M，M表示样品在360度旋转范围内的投影总数目。

所述

步骤三、获取梯度角及梯度角直方图；

求取合成投影P^ave的所有像素点的梯度角，P^ave内任意投影值p(x_i,y_j)梯度角β(x_i,y_j)的计算公式如下：

其中，表示任意投影值p(x_i,y_j)在X_d轴方向的梯度值，表示任意投影值p(x_i,y_j)在Y_d轴方向的梯度值。β表示任意像素的梯度角，x_i为探测器坐标系X_d轴上的变量，y_j为探测器坐标系Y_d轴上的变量。

统计所有像素点的梯度角的直方图，梯度角直方图的计算式如下：

其中F(β)表示直方图函数，该函数自变量为像素点的梯度角β，w表示合成投影图像P^ave的宽度，h表示合成投影图像P^ave的高度，n(β)表示梯度角等于β的像素点的个数。

步骤四、获取旋转轴线偏摆角；

求取梯度角直方图函数F(β)的中心点坐标β₀，该坐标即为旋转轴投影线41与X_d轴的夹角α，计算公式如下：

步骤五、旋转轴线偏摆角误差补偿；

将计算出的旋转轴投影线41与X_d轴的夹角α作为修正参数引入到中子层析成像系统的重建算法中，即可实现对旋转轴线偏摆角的误差补偿。

2.根据权利要求1所述的一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法，其特征在于：能够用于补偿旋转轴投影线41与X_d轴的夹角不等于90度的旋转轴线偏摆角的误差补偿。

3.根据权利要求1或2所述的一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法，其特征在于：被检测样品3在旋转平台2的带动下旋转360度，旋转平台2步进旋转，步进角步长一般取0.5°～1°。

4.根据权利要求1或2所述的一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法，其特征在于：投影总数目取值为720～1200。

5.根据权利要求1或2所述的一种适用于中子层析成像系统的旋转轴线偏摆角的误差补偿方法，其特征在于：合成投影图像P^ave中的像素点A与投影像素点B之间存在的关系为：

。