用于CT系统中标定几何参数的标定体模、标定装置及标定方法

摘要

本发明涉及一种用于CT系统中几何参数标定的标定体模，包括两个相同的平行设置的标定板及小球，所述两个标定板上分别对应开设有若干相同的小孔，所述小球至少为四个且分别安装于所述两个标定板上的小孔内。本发明涉及一种用于CT系统中几何参数标定的标定装置包括由光源、旋转台和探测器组成的CT系统，所述标定装置还包括垂直放置在所述光源和探测器之间的旋转台上的标定体模。另外，本发明还提供了一种用于CT系统中几何参数标定的标定方法。上述标定装置及标定方法仅需对标定体模进行1次投影就能够同时标定CT系统的全部7个几何参数，操作简单，极大满足了后续CT系统的图像校正。

说明书

技术领域

本发明涉及CT系统中几何参数的标定，尤其涉及一种用于CT系统中几何参数标定的标定体模、标定装置及标定方法。

背景技术

计算机断层成像(CT)是通过无损方式获取物体内部结构信息的一种重要成像手段，它拥有高分辨率、高灵敏度以及多层次等众多优点，被广泛应用于各个医疗临床检查领域。几何参数标定是显微CT系统调试的重要组成部分，也是获得精确重建CT图像的前提条件。

中国专利申请号201210148432.8提出的“一种X射线锥束计算机层析成像系统的几何参数标定方法”。该方法通过相机标定技术来复现锥束CT系统中的X射线源、平板探测器和转轴之间的几何位置关系，从而对锥束CT的几何参数及其误差进行直接求解，是一种系统化的测量求解方法，可以将锥束CT抽象为基本的相机系统模型，从而同时求解系统中相互关联的多个几何参数。

但是，上述公开的CT系统的几何参数标定方法不能同时标定CT系统的全部7个几何参数；且并不能同时适用于离体CT系统(标定体模旋转，光源和探测器不动)和活体CT系统(标定体模固定，光源和探测器旋转)；另外，参数标定过程较复杂，需要对标定体模进行多次投影，难于多次重复实现。

发明内容

基于此，有必要针对上述CT系统中几何参数的标定系统存在的缺陷，提供一种能够同时标定CT系统中几何参数的标定体模。

一种用于CT系统中几何参数标定的标定体模，包括：包括两个相同的平行设置的标定板及小球，所述两个标定板上分别对应开设有若干相同的小孔，所述小球至少为四个且分别安装于所述两个标定板上的小孔内。

本发明还提供了一种用于CT系统中几何参数标定的标定装置，包括由光源、旋转台和探测器组成的CT系统，所述标定装置还包括垂直放置在所述光源和探测器之间的旋转台上的标定体模，所述标定体模包括两个相同的平行设置的标定板及小球，所述两个标定板上分别对应开设有若干相同的小孔，所述小球至少为四个且分别安装于所述两个标定板上的小孔内；所述标定体模随所述旋转台旋转，所述探测器用于采集各个小球的投影图像的中心点坐标；及图像处理模块，用于根据所述投影图像的中心点坐标构建第一投影矩阵，所述图像处理模块还用于根据所述CT系统的几何结构及所述标定体模几何结构分别计算各个小球在所述CT系统坐标系中的空间位置坐标，以及根据所述空间位置坐标构建第二投影矩阵，根据所述第一投影矩阵及所述第二投影矩阵计算得到探测器平面的投影中心点坐标、光源到探测器的距离、探测器的扭转角、探测器的倾斜角、探测器的旋转角及光源到旋转台旋转轴的距离。

另外，本发明还提供了一种CT系统中几何参数的标定方法，所述CT系统包括光源、旋转台和探测器，所述标定方法包括下述步骤：在所述光源和探测器之间的旋转台上垂直放置标定体模，所述标定体模包括两个相同的平行设置的标定板及小球，所述两个标定板上分别对应开设有若干相同的小孔，所述小球至少为四个且分别安装于所述两个标定板上的小孔内；将所述标定体模随所述旋转台旋转，通过所述探测器采集各个小球的投影图像的中心点的横坐标和纵坐标；根据所述投影图像的中心点坐标构建第一投影矩阵；根据所述CT系统的几何结构关系和标定体模的几何结构关系分别计算各个小球在所述CT系统坐标系中的空间位置坐标；根据所述空间位置坐标构建第二投影矩阵；根据所述第一投影矩阵及所述第二投影矩阵计算得到探测器平面的投影中心点坐标、光源到探测器的距离、探测器的扭转角、探测器的倾斜角、探测器的旋转角及光源到旋转台旋转轴的距离。

上述用于CT系统中几何参数的标定装置及标定方法，采用简单的标定体模，在旋转台上旋转，通过探测器采集标定体模上各个小球的投影图像的中心点坐标，根据投影图像的中心点坐标构建第一投影矩阵，根据CT系统的几何结构及标定体模几何结构分别计算各个小球在所述标定装置坐标系中的空间位置坐标，并根据所述空间位置坐标构建第二投影矩阵，根据第一投影矩阵及第二投影矩阵计算得到探测器平面的投影中心点坐标、光源到探测器的距离、探测器的扭转角、探测器的倾斜角、探测器的旋转角及光源到旋转台旋转轴的距离。

上述标定装置及标定方法仅需对标定体模进行1次投影就能够同时标定CT系统的全部7个几何参数，操作简单，极大满足了后续CT系统的图像校正；另外，上述标定装置及标定方法可以同时适用于离体CT系统和活体CT系统，适应性广；同时，上述标定体模结构简单、制备成本低廉，易于获得。

附图说明

图1为本发明实施例一提供的用于CT系统中几何参数标定的标定体模的示意图。

图2为本发明提供的实施例一提供的标定板上开设的小孔的结构示意图。

图3为本发明实施例二提供的用于CT系统中几何参数标定的标定装置的结构示意图。

图4为本发明实施例三提供的CT系统中几何参数的标定方法的流程图。

图5为本发明提供的实施例三提供的通过探测器采集各个小球的投影图像的中心点坐标的步骤流程图。

图6为本发明实施例三提供的根据CT系统的几何结构关系和标定体模的几何结构关系分别计算各个小球在CT系统坐标系中的空间位置坐标步骤流程图。

图7为本发明实施例四提供的当标定体模上嵌入18个金属小球时得到的投影图像。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及具体实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例一：

请参阅图1及图2，为本发明实施例一提供的用于CT系统中几何参数标定的标定体模的示意图。

标定体模100包括标定板110及小球120。

标定板110的数量为两个且结构相同。两个标定板110正对且平行设置。在本发明提供的实施例中，两个标定板的间距优选为46mm，每个标定板的宽度优选为150mm、厚度优选为2mm。两个标定板110上分别对应开设有若干相同的小孔111。小孔111呈方阵式分布于标定板上。请参阅图2，为本发明提供的实施例一提供的标定板110上开设的小孔111的结构示意图，小孔111呈7*7的方阵分布于标定板110上，且小孔111的直径优选为1mm、孔深度优选为0.5mm，相邻孔之间间距优选为20mm。

小球120至少为四个。小球120分别安装于两个标定板110上的小孔111内，即任意一个标定板110上均安装有小球120。可以理解，小孔120的数量还可以为6个或其他更多的数量。

在本发明提供的实施例中，标定板110的材质优选为聚乙烯塑料。小球120优选为不锈钢小球。

实施例二：

请参阅图3，图3为本发明实施例二提供的用于CT系统中几何参数标定的标定装置的结构示意图。

标定装置200包括由光源210、旋转台(图未示)和探测器230组成的CT系统、标定体模100及图像处理模块(图未示)。

光源210和探测器230之间设有旋转台。标定体模100垂直设置于旋转台上。标定体模100随旋转台旋转。探测器230用于采集标定体模100上开设的小孔111内的各个小球120的投影图像的中心点坐标。可以理解，小球120在成像平面上的投影图像的中心点坐标均不重合。

图像处理模块用于根据投影图像的中心点坐标构建第一投影矩阵。图像处理模块还用于根据CT系统的几何结构及标定体模100几何结构分别计算各个小球120在CT系统坐标系中的空间位置坐标，以及根据空间位置坐标构建第二投影矩阵。图像处理模块根据第一投影矩阵及第二投影矩阵计算得到探测器230平面的投影中心点坐标、光源210到探测器230的距离、探测器230的扭转角、探测器230的倾斜角、探测器230的旋转角及光源210到旋转台旋转轴的距离。

实施例三：

请一并参阅图3及图4，其中图4为本发明实施例三提供的CT系统中几何参数的标定方法的流程图，其中CT系统包括光源210、旋转台和探测器230，具体包括下述步骤：

步骤S310：在光源210和探测器230之间的旋转台上垂直放置标定体模100。在本发明提供的实施例中，标定体模100与本发明实施例一提供的标定体模相同。可以理解，小球120至少为四个且分别安装于两个标定板110上的小孔111内，即标定体模100的两个标定板110上均要放置小球120。

步骤S320：将标定体模100随旋转台旋转，通过探测器230采集各个小球120的投影图像的中心点的横坐标和纵坐标。

请参阅图5，为本发明提供的实施例三提供的通过探测器230采集各个小球的投影图像的中心点坐标的步骤流程图，其包括下述步骤：

步骤S321：定义探测器230所在的平面作为坐标平面。其中以探测器230左下角为坐标原点，探测器230上像素的排列方向U、V分别为横坐标轴和纵坐标轴。

步骤S322：将标定体模100随旋转台旋转。可以理解，标定体模100上开设的小孔111内的小球120随旋转台匀速旋转，通过光源210对小球120进行投影，并得到投影图像。在本发明提供的实施例中，投影的次数优选为1次，可以理解，投影的次数还可以为多次，当投影的次数越多时，最终得到的CT系统的几何参数越准确。

步骤S323：探测器230采集各个小球120的投影图像。可以理解，每次投影小球都在探测器230上获得1个投影图像，即小球都在以探测器230所在坐标平面为成像平面上得到投影图像。

步骤S324：通过边缘提取处理方法依次获得各个小球120投影图像的中心点在以探测器230所在的平面作为坐标平面的横坐标和纵坐标。可以理解，小球120在成像平面上的投影图像的中心点坐标均不重合。

步骤S330：根据投影图像的中心点横坐标和纵坐标构建第一投影矩阵。在本发明提供的实施例中，根据投影图像的中心点横坐标和纵坐标构建第一投影矩阵具体如下：

将各个小球120的投影图像的中心点的横坐标和纵坐标分别代入矩阵[u，v，1]^T中，构建第一投影矩阵[u，v，1]^T，其中，u和v分别表示小球120投影中心点的横坐标和纵坐标，[u，v，1]^T为3xN的矩阵，N为小球个数。

步骤S340：根据CT系统的几何结构关系和标定体模的几何结构关系分别计算各个小球在CT系统坐标系中的空间位置坐标。

请参阅图6，为本发明实施例三提供的根据CT系统的几何结构关系和标定体模100的几何结构关系分别计算各个小球120在CT系统坐标系中的空间位置坐标步骤流程图，具体包括下述步骤：

步骤S341：定义坐标系(X，Y，Z)，其中Y轴为旋转台的旋转轴，Z轴为光源210到旋转轴的射线。

步骤S342：使标定体模100的中心点所在的轴线与旋转轴重合。在本发明提供的标定体模100的另一实施例中，沿标定体模100的中心点所在的轴线上设有一金属圆柱，该金属圆柱与旋转轴重合，且该金属圆柱能绕旋转轴旋转。

步骤S343：旋转旋转台使两个标定板110与探测器230所在的平面平行，即可获取各个小球在坐标系(X，Y，Z)上的空间位置坐标(x，y，z)。可以理解，由于标定体模100的结构是已知的，这样标定体模100上的各个小球120的相对空间位置确定，将标定体模100固定在旋转台上，并使标定体模100中金属圆柱与旋转轴Y轴重合。旋转旋转台，使得标定体模100的两个标定板110所在的平面与探测器230所在的平面平行，光源210在标定体模100的两个标定板110的垂点分别为两个平面的中心，即可得到不锈钢小球120在坐标系(X，Y，Z)上的空间位置坐标(x，y，z)。

步骤S350：根据空间位置坐标构建第二投影矩阵。在本发明提供的实施例中，其中，根据不锈钢小球120在坐标系(X，Y，Z)上的空间位置坐标(x，y，z)构建第二投影矩阵，具体包括下述步骤：将各个小球的空间位置坐标(x，y，z)分别代入矩阵[x，y，z，1]^T中，构建第二投影矩阵[x，y，z，1]^T。其中，x，y，z为各个小球的空间位置坐标(x，y，z)。

步骤S360：根据第一投影矩阵及第二投影矩阵计算得到探测器平面的投影中心点坐标、光源到探测器的距离、探测器的扭转角、探测器的倾斜角、探测器的旋转角及光源到旋转台旋转轴的距离。

根据第一投影矩阵及第二投影矩阵计算得到探测器平面的投影中心点坐标、光源到探测器的距离、探测器的扭转角、探测器的倾斜角、探测器的旋转角及光源到旋转台旋转轴的距离的步骤流程图，具体包括下述步骤：

定义坐标系(X，Y，Z)，其中Y轴为旋转台的旋转轴，Z轴为光源210到旋转轴的射线。

定义探测器230左下角为坐标原点，U、V为探测器230上像素的排列方向，λ为探测器像素的大小、(u₀、v₀)为探测器230平面的投影中心点坐标，D为光源210到探测器230的距离，α为探测器230的扭转角、β为探测器230的倾斜角、γ为探测器230的旋转角，R为光源210到旋转轴的距离。

根据第一投影矩阵[u，v，1]^T及第二投影矩阵[x，y，z，1]^T构建公式a[u，v，1]^T＝A[x，y，z，1]^T，其中，A为系统矩阵，a为系数。可以理解，a可以根据经验选取。

根据公式A＝Φ[ζ|ω]并计算Φ、ζ，其中，

$>\phi =\left(\begin{array}{ccc}\frac{D}{\lambda }& 0& u_0\\ 0& \frac{D}{\lambda }& v_0\\ 0& 0& 1\\ & & \end{array}\right)$>

$>\zeta =\left(\begin{array}{ccc}\cos \beta \cos \gamma & \sin \alpha \sin \beta \cos \gamma -\cos \alpha \sin \gamma & \sin \alpha \sin \gamma +\cos \alpha \sin \beta \cos \gamma \\ \cos \beta \sin \gamma & \cos \alpha \cos \gamma +\sin \alpha \sin \beta \sin \gamma & \cos \alpha \sin \beta \sin \gamma -\sin \alpha \cos \gamma \\ -\sin \beta & \sin \alpha \cos \beta & \cos \alpha \cos \beta \end{array}\right)$>

$>\omega =\left(\begin{array}{c}{\omega }_1\\ {\omega }_2\\ {\omega }_3\end{array}\right)$>

ω₁＝A₃₄、ω₂＝(A₂₄-Φ₂₃A₃₄)/Φ₂₃、ω₂＝(A₂₄-Φ₂₃A₃₄)/Φ₂₃、ω₂＝(A₁₄-Φ₁₃A₃₄-ω₂A₁₂)/Φ₁₁。上述公式中，D为光源到探测器的距离，(u₀、v₀)为光源210在探测器230平面的投影中心点坐标，λ为探测器230像素的大小，α为探测器230的扭转角、β为探测器230的倾斜角、γ为探测器230的旋转角。

可以理解，先获得N个小球120在坐标系(X，Y，Z)上的空间位置坐标(x_i，y_i，z_i)，i＝1，2，3，......N和以探测器230所在的平面作为坐标平面的横坐标和纵坐标(u_i，v_i)，i＝1，2，3，......N。由公式a[u，v，1]^T＝A[x，y，z，1]^T，有如下方程式：

u_ia_i＝A₁₁x_i+A₁₂y_i+A₁₃z_i+A₁₄

v_ia_i＝A₂₁x_i+A₂₂y_i+A₂₃z_i+A₂₄

ω_i＝A₃₁x_i+A₃₂y_i+A₃₃z_i+A₃₄

消掉系数a_i，得到如下两个方程：

A₁₁x_i+A₁₂y_i+A₁₃z_i+A₁₄-u_i(A₃₁x_i+A₃₂y_i+A₃₃z_i+A₃₄)＝0

A₂₁x_i+A₂₂y_i+A₂₃z_i+A₂₄-v_i(A₃₁x_i+A₃₂y_i+A₃₃z_i+A₃₄)＝0

将以上方程写成矩阵形式，得到：H·A*＝0，其中，

$>H=\left(\begin{array}{cccccccccccc}{x}_i& y_i& z_i& 1& 0& 0& 0& 0& -u_i{x}_i& {-u}_i{y}_i& -u_i{z}_i& -u_i\\ 0& 0& 0& 0& x_i& y_i& z_i& 1& -v_i{x}_i& -v_i{y}_i& -v_i{z}_i& -v_i\end{array}\right)$>

$>A^{*}=\left(\begin{array}{c}{A}_{11}\\ A_{12}\\ A_{13}\\ A_{14}\\ A_{21}\\ A_{22}\\ A_{23}\\ A_{24}\\ A_{31}\\ A_{32}\\ A_{33}\\ A_{34}\end{array}\right)$>

对于方程H·A*＝0，可以通过奇异值分解方法，对矩阵H进行分解，求出A*，即可得系统矩阵A。

根据公式A＝Φ[ζ|ω]，可以得到A_3×3＝Φζ(即去掉了矩阵A的3x4的最后一列)，由于Φ为上三角矩阵，可以对矩阵A_3×3进行“QR分解”来得到Φ，ζ。

依据公式o＝[o_x，o_y，o_z]^T＝-ζ^Tω，计算XYZ坐标系中心点坐标o_x、o_y及o_z。

其中，上述式中，ζ可由上述步骤求出，

$>\omega =\left(\begin{array}{c}{\omega }_1\\ {\omega }_2\\ {\omega }_3\end{array}\right)$>

ω₁＝A₃₄、ω₂＝(A₂₄-Φ₂₃A₃₄)/Φ₂₃、ω₂＝(A₂₄-Φ₂₃A₃₄)/Φ₂₃、ω₂＝(A₁₄-Φ₁₃A₃₄-ω₂A₁₂)/Φ₁₁。

计算探测器230平面的投影中心点坐标(u₀、v₀)，光源210到探测器230的距离D，探测器230的扭转角α、探测器230的倾斜角β、探测器230的旋转角γ，光源210到旋转轴的距离R。

可以理解，依据上述得到的Φ，ζ及XYZ坐标系中心点坐标o_x、o_y及o_z，可以计算得到探测器230平面的投影中心点坐标(u₀、v₀)，光源210到探测器230的距离D，探测器230的扭转角α、探测器230的倾斜角β、探测器230的旋转角γ，光源210到旋转轴的距离R。

其中，探测器230平面的投影中心点坐标(u₀、v₀)的计算公式为：u₀＝Φ₁₃，v₀＝Φ₂₃，Φ₁₃表示矩阵Φ的第1行和第3列对应的元素，Φ₂₃表示矩阵Φ的第1行和第3列对应的元素。

光源210到探测器230的距离D的计算公式为：D＝Φ₁₁λ，Φ₁₁表示矩阵Φ的第1行和第1列对应的元素。

探测器230的扭转角α的计算公式为：α＝Arctan2(ζ₃₂，ζ₃₃)，ζ₃₂表示矩阵ζ的第3行和第2列对应的元素，ζ₃₃表示矩阵ζ的第3行和第3列对应的元素。

探测器230的倾斜角β的计算公式为：β＝sin^-1(-ζ₃₁)，ζ₃₁表示矩阵ζ的第3行和第1列对应的元素。

探测器230的旋转角γ的计算公式为：γ＝Arctan2(ζ₁₂，ζ₃₁)，ζ_，12表示矩阵ζ的第1行和第2列对应的元素，ζ₃₁表示矩阵ζ的第，3行和第1列对应的元素。

光源210到旋转轴的距离R的计算公式为：o_x、o_y、o_z为XYZ坐标系中心点坐标。

实施例四：

请参阅图7，为本发明实施例四提供的当标定体模上嵌入18个金属小球时得到的投影图像。在标定体模100上嵌入18个金属小球时的情况。分别在上下两块标定板上各对应放入9个小球，得到的投影图像。为了验证算法的有效性，按极端情况，只进行了1次投影(等于1)。理论上，当取值越大时，结果会越接近于真实值。依据上述的标定方法计算仿真得到的系统几何参数模拟值与真实值如表格1所示。从结果可以看出，该方法精度非常高。

表格1

RDαβγu₀v₀真实值5701005-1.5708-1.57080.00257257计算值571.77601008.5622-1.5708-1.57050.00257255.0620误差0.312％0.354％0％0.019％0％0％0.754％

上述CT系统中几何参数的标定方法和标定装置采用标定体模旋转，光源和探测器不动，从而实现CT系统中几何参数的标定；可以理解，对于活体CT系统(标定体模固定，光源和探测器旋转)同样可以实现CT系统中几何参数的标定。

上述用于CT系统中几何参数的标定装置及标定方法，采用简单的标定体模，在旋转台上旋转，通过探测器采集标定体模上各个小球的投影图像的中心点坐标，根据投影图像的中心点坐标构建第一投影矩阵，根据CT系统的几何结构及标定体模几何结构分别计算各个小球在所述标定装置坐标系中的空间位置坐标，并根据所述空间位置坐标构建第二投影矩阵，根据第一投影矩阵及第二投影矩阵计算得到探测器平面的投影中心点坐标、光源到探测器的距离、探测器的扭转角、探测器的倾斜角、探测器的旋转角及光源到旋转台旋转轴的距离。

上述标定装置及标定方法仅需对标定体模进行1次投影就能够同时标定CT系统的全部7个几何参数，操作简单，极大满足了后续CT系统的图像校正；另外，上述标定装置及标定方法可以同时适用于离体CT系统和活体CT系统，适应性广；同时，上述标定体模结构简单、制备成本低廉，易于获得。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。