CN102103757B - 锥束图像重建方法及装置 - Google Patents

Abstract

一种锥束图像重建方法，包括如下步骤：扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵；根据所述三维运动矩阵，对非等中心投影图像转化，得到等中心投影图像；根据所述等中心投影图像，重建锥束图像。上述锥束图像重建方法及装置通过跟踪C型臂的运动姿态，将扫描得到的非等中心投影图像结合投影的几何关系，在不增加额外跟踪设备的基础上，转化为等中心图像，有效地降低了成本，提高投影图像的准确性，提高了投影图像的成像质量及其鲁棒性。

Description

锥束图像重建方法及装置

【技术领域】

本发明涉及计算机成像领域，特别是涉及一种锥束图像重建方法及装置。

【背景技术】

X射线断层扫描计算机成像（computed tomography，简称CT）是计算机层析成像技术，不仅应用于放射诊断医学领域，同时也是现代工业无损检测和勘探领域的重要手段，它是利用射线束对物体进行扫描，根据扫描所得结果运用特定的重建算法，以二维或三维成像的方式呈现出物体内部的密度分布。

为此人们提出了许多图像重建的方法，大致可以分为两大类，即解析法和迭代法。解析法在实际的应用中实现的条件过于严格，后者如联合代数重建算法（SART）、最大似然-期望算法（ML-EM）以及有序子集-期望值最大化算法（OS-EM），与解析法相比较，迭代法更适合应用于有限角投影和噪声状况较差等情况，但是在实际的应用过程中，又存在着重建速度慢、计算量大、鲁棒性不好及重建效率不高的问题。

此外，通过传统的C型臂所获得的x射线图片由于C型臂旋转中心的偏移，导致x射线图片是非等中心的锥束图像，若直接用非等中心图像进行重建，重建结果比较差，因此有必要对其非等中心进行转化。

【发明内容】

基于此，有必要提供一种提高图像质量的锥束图像重建方法。

此外，还有必要提供一种提高图像质量的锥束图像重建装置。

一种锥束图像重建方法，包括如下步骤：扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵；所述扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵的具体过程为：

提取体模基准点，根据所述体模基准点进行C型臂的运动姿态估计，运用投影图像上至少3组基准点建立方程组得到三维运动矩阵；

所述姿态估计的步骤是：

选取基准点，设定三维运动矩阵初值，建立残差方程计算得到所述基准点的初始旋转分量和初始平移分量，残差方程如下所示：

$\mathrm{\δ}{x}_i=x_i(R_k,T_k)-{\stackrel{\‾}{x}}_i$

$\mathrm{\δ}{y}_i=y_i(R_k,T_k)-{\stackrel{\‾}{y}}_i$

其中，(x_i,y_i)为基准点的投影坐标，

为基准点的标注，(δx_i,δy_i)为残差值，R_k为F_x的旋转矩阵，T_k为F_x的平移向量，F_x为三维运动矩阵；

再次选取基准点，将初始的旋转分量和平移分量进行迭代，得到旋转增量及平移增量，并更新得到旋转分量和平移分量；

判断相邻的两次选取基准点所得到的旋转分量之差及平移分量之差是否小于预设阈值或所述迭代次数是否大于预设次数，是，则终止迭代，反之，则再次选取基准点进行迭代；根据所述三维运动矩阵，对非等中心投影图像转化，得到等中心投影图像；根据所述等中心投影图像，重建锥束图像。

优选地，所述提取体模基准点的步骤是：对所述投影图像进行高斯滤波；将滤波后得到的投影图像进行二值化处理；通过连通成分标记进行基准点的标注。

优选地，所述提取体模基准点的步骤是：剔除投影图像中CT值大于第一阈值的区域，得到初选的投影区域；在所述初选的投影区域中提取CT值大于第二阈值的感兴趣的投影区域，所述第一阈值大于第二阈值；通过连通成分标记进行基准点的标注。

优选地，所述迭代次数至少1次。

优选地，所述根据所述三维运动矩阵，对非等中心投影图像转化，得到等中心投影图像的步骤是：通过投影关系得到投影图像的投影坐标，根据旋转分量和平移分量计算得到圆心偏移量，并根据所述圆心偏移量和锥束射线的扫描角度，还原所述投影图像的物理坐标；转换所述物理坐标得到等中心的投影图像。

优选地，所述根据所述等中心投影图像，重建锥束图像的步骤是：通过等中心投影图像中体素的衰减系数和所述投影图像获得体素值向量；修正所述等中心投影图像，得到修正值；根据所述修正值和体素的衰减系数进行反投影，更新体素值向量。

一种锥束图像重建装置，至少包括：跟踪模块，用于扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵；所述跟踪模块还用于提取体模基准点，根据所述体模基准点进行C型臂的运动姿态估计，运用投影图像上至少3组基准点建立方程组得到三维运动矩阵；

进一步的，所述跟踪模块选取基准点，设定三维运动矩阵初值，建立残差方程计算得到所述基准点的初始旋转分量和初始平移分量，残差方程如下所示：

$\mathrm{\δ}{x}_i=x_i(R_k,T_k)-{\stackrel{\‾}{x}}_i$

$\mathrm{\δ}{y}_i=y_i(R_k,T_k)-{\stackrel{\‾}{y}}_i$

其中，(x_i,y_i)为基准点的投影坐标，

为基准点的标注，(δx_i,δy_i)为残差值，R_k为F_x的旋转矩阵，T_k为F_x的平移向量，F_x为三维运动矩阵；再次选取基准点，将初始的旋转分量和平移分量进行迭代，得到旋转增量及平移增量，并更新得到旋转分量和平移分量；判断相邻的两次选取基准点所得到的旋转分量之差及平移分量之差是否小于预设阈值或所述迭代次数是否大于预设次数，是，则终止迭代，反之，则再次选取基准点进行迭代；中心变换模块，用于根据所述三维运动矩阵，对非等中心投影图像转化，得到等中心投影的投影图像；重建模块，用于根据所述等中心投影图像，重建锥束图像。

优选地，所述提取单元对所述投影图像在高斯滤波后进行二值化处理，并通过连通成分标记进行基准点的标注。

优选地，所述提取单元剔除投影图像中CT值大于第一阈值的区域，得到初选的投影区域，在所述初选的投影区域中提取CT值大于第二阈值的感兴趣的投影区域，所述第一阈值大于第二阈值，并通过连通成分标记进行基准点的标注。

优选地，所述估计单元中的迭代次数至少1次。

优选地，所述中心变换模块包括：还原单元，用于通过投影关系得到非等中心投影图像的投影坐标，根据旋转分量和平移分量计算得到圆心偏移量。并根据圆心偏移量和锥束射线的扫描角度，还原投影图像的物理坐标；转换单元，用于转换所述物理坐标得到等中心的投影图像。

优选地，所述重建模块包括：投影单元，用于通过等中心投影图像中体素的衰减系数和所述投影图像获得体素值向量；修正单元，用于修正所述等中心投影图像，得到修正值；反投影单元，用于根据所述修正值和体素的衰减系数进行反投影，更新体素值向量。

上述锥束图像重建方法及装置通过跟踪C型臂的运动姿态，将扫描得到的非等中心投影图像结合投影的几何关系，在不增加额外跟踪设备的基础上，转化为等中心图像，有效地降低了成本，提高投影图像的准确性，提高了投影图像的成像质量及其鲁棒性。

上述锥束图像重建方法及装置通过多次运动姿态估计，来得到三维运动矩阵，收敛性好，有效地降低了投影点与基准点坐标之间的误差。

【附图说明】

图1为一实施例中锥束图像重建方法的流程图；

图2为一实施例中步骤S10的方法流程图；

图3为一实施例中体模的示意图；

图4为一实施例中提取体模基准点的方法流程图；

图5为另一实施例中提取体模基准点的方法流程图；

图6为一实施例中姿态估计的方法流程图；

图7为一实施例中步骤S20的方法流程图；

图8为一实施例中非等中心投影的几何关系图；

图9为一实施例中等中心投影的几何关系图；

图10为一实施例中重建锥束图像的方法流程图；

图11为一实施例中锥束图像重建装置的模块图。

图12为一实施例中原始等中心投影图像与等中心图像差；

图13为一实施例中等中心投影图像示例；

图14为一实施例中非等中心投影图像示例；

图15a、图15b和图15c为一实施例中脊骨原始三维截面图；

图16a、图16b和图16c为一实施例中重建后脊骨三维截面图。

【具体实施方式】

图1示出了一实施例中的锥束图像重建的方法流程，包括如下步骤：

在步骤S10中，扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵。一实施方式中，C型臂中所扫描得到的投影图像是非等中心的，即投影图像的中心与实际的旋转中心存在着偏差，因此需要跟踪C型臂的运动姿态以修正非等中心投影。如图2所示，步骤S10的过程具体是：

在步骤S 120中，提取体模基准点。一实施例中，如图3所示，体模为镶嵌预设数量铅球的物体，该铅球呈螺旋形排列，因此在其投影上表现为若干像素的阴影。体模中的铅球质地均匀，内部没有瑕疵，铅球半径在1.5~3mm内，间距为3~8mm，从而保证铅球的投影不会产生重叠。具体地，如图4所示，提取体模基准点的过程是：

在步骤S122中，对投影图像进行高斯滤波。为排除噪声对投影图像的影响，通过高斯滤波器对投影图像进行预处理。

在步骤S124中，将滤波后得到的投影图像进行二值化处理。二值化处理是把图像分成前景和背景的方法，基于滤波后所得到的投影图像具有较高对比度的情况，可以通过图像二值化的方法来进行基准点的分割。

假定原灰度图像为f(x,y)，则二值化后的新图像为：

$p(x,y)=\left\{\begin{array}{c}1,f(x,y)\&GreaterEqual;T\\ 0,f(x,y)<T\end{array}\right.$

其中，T为f(x,y)k领域的平均值，

其中，k为阈值，M为k邻域内像素个数。

在步骤S126中，通过连通成分标记进行基准点的标注。连通成分标记是逐个进行体素扫描，扫描方式是沿图像列顺序遍历，然后沿图像行顺序遍历。查找体素的相邻体素，如果体素不是背景且没有相邻体素，则标识该体素，并继续遍历，如果体素不是背景并且有相邻的体素，该体素和相邻的体素都分配了不同的标记数，则将该体素与其相邻的体素都标示为同一连通成分，即把该体素与其相邻体素的所有标记数中的最小值分配给这一连通成分，并存储其相邻关系及标识数，最后用最小的标识数重新标识体素。通过连通成分标记既可以得到所有的连通成分，又可以在此基础上独立计算连通区域的统计特性。由于体模中的基准点在投影图像上为若干像素的集合，形心为体素的几何中心，因此对每一基准点均用像素的形心(i,j)来表示，具体公式如下所示：

$\stackrel{\&OverBar;}{i}=\frac{\underset{i=r}{\overset{m-1}{\mathrm{\&Sigma;}}}i\&times;f(x,y)}{\underset{i=r}{\overset{m-1}{\mathrm{\&Sigma;}}}f(x,y)}$ $\stackrel{\&OverBar;}{j}=\frac{\underset{j=c}{\overset{n-1}{\mathrm{\&Sigma;}}}j\&times;f(x,y)}{\underset{j=c}{\overset{n-1}{\mathrm{\&Sigma;}}}f(x,y)}$

其中i，j所选用的区域为基准点所在的联通区域。

其他实施例中，体模的投影表现为若干像素、高对比度的投影，因此也可以通过阈值的方法提取体模基准点。为增加所提取到的基准点的鲁棒性，采用双层提取的方式，如图5所示，提取体模基准点的过程是：

在步骤S132中，剔除投影图像中CT值大于第一阈值的区域，得到初选的投影区域。具体地，第一阈值优选为2800~3200Hounsfield uints。Hounsfield uint为CT值的单位，简称HU。CT值是代表射线被吸收后的衰减值，人体组织、金属等对射线的吸收各有不同，金属的衰减值较大，则在投影图像上更暗，因此可以通过设定的第一阈值来提取铅球的投影，基于铅球的尺寸、形状及体模上铅球的排列位置，采用尺寸和距离约束，以保证提取出的投影为铅球投影。

在步骤S134中，在被选的投影区域中提取CT值大于第二阈值的感兴趣的投影区域，且第一阈值大于第二阈值。具体地，第二阈值优选为1800~2200Hounsfield uints。在初选得到的投影区域中用一较小的阈值，即CT值大于2000Hounsfield uints时，提取感兴趣的投影区域，同样，与第一次提取相同，也对所提取的铅球投影进行尺寸和距离的约束，从而保证投影结果。

在步骤S136中，通过连通成分标记进行基准点标注。如前所述，连通成分标记是逐个体素扫描图像，扫描方式是沿图像列顺序遍历，然后沿图像行顺序遍历。查找体素的相邻体素，如果体素不是背景且没有相邻体素，则标识该体素，并继续遍历，如果体素不是背景并且有相邻的体素，该体素和相邻的体素都分配了不同的标记数，则将该体素与其相邻的体素都标示为同一连通成分，即把该体素与其相邻体素的所有标记数中的最小值分配给这一连通成分，并存储其相邻关系及标识数，最后用最小的标识数重新标识体素。通过连通成分标记既可以得到所有的连通成分，又可以在此基础上独立计算连通区域的统计特性。由于体模中的基准点在投影图像上为若干像素的集合，因此对每一基准点均用像素的形心来表示，具体公式如下所示：

$\stackrel{\&OverBar;}{i}=\frac{\underset{i=r}{\overset{m-1}{\mathrm{\&Sigma;}}}i\&times;f(x,y)}{\underset{i=r}{\overset{m-1}{\mathrm{\&Sigma;}}}f(x,y)}$ $\stackrel{\&OverBar;}{j}=\frac{\underset{j=c}{\overset{n-1}{\mathrm{\&Sigma;}}}j\&times;f(x,y)}{\underset{j=c}{\overset{n-1}{\mathrm{\&Sigma;}}}f(x,y)}$

在步骤S140中，根据体模基准点，进行C型臂的运动姿态估计，得到三维运动矩阵。由于旋转中心偏移，使得投影图像是非等中心的，因此需要跟踪C型臂的空间三维情况。C型臂的三维运动可以由六个自由量所确定，即平移分量和旋转分量。由射线方程组可以知道，为求出这六个自由量，理想的情况是运用投影图像上三组对应基准点来建立方程组，但是考虑到投影图像的变形以及鲁棒性，所选择的基准点数量为至少3组。

一实施方式中，如图6所示，姿态估计的步骤具体是：

在步骤S142中，选取基准点，设定三维运动矩阵初值，建立残差方程计算得到基准点初始旋转分量和初始平移分量。一实施例中，为减少计算处理量，加快图像重建速度，间隔选取基准点，在优选的实施例中，体模中设定了15个基准点，每间隔两个基准点选取一个，共选取5个基准点。三维运动矩阵为

$F_x=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& r_{14}\\ r_{21}& r_{22}& r_{23}& r_{24}\\ r_{31}& r_{32}& r_{33}& r_{34}\\ r_{41}& r_{42}& r_{43}& r_{44}\end{array}\right]$

其中，r_ij为旋转分量，T_ij为平移分量，F_x为变换矩阵，表示物体相对于射线源的变换，三维运动矩阵的假定初值为F_x＝[1000 0100 0010 0001]。

建立残差方程，并将基准点的投影坐标及标注导入残差方程中，得到初始旋转分量和初始平移分量，残差方程如下所示：

$\mathrm{\&delta;}{x}_i=x_i(R_k,T_k)-{\stackrel{\&OverBar;}{x}}_i$

$\mathrm{\&delta;}{y}_i=y_i(R_k,T_k)-{\stackrel{\&OverBar;}{y}}_i$

其中，(x_i,y_i)为基准点的投影坐标，

为基准点的标注，(δx_i,δy_i)为残差值，R_k为F_x的旋转矩阵，T_k为F_x的平移向量。

在残差方程中，根据体模中基准点的投影坐标及与标注之间的关系，反求出当前的物体相对于射线源的三维变化。

在步骤S144中，再次选取基准点，将初始的旋转分量和平移分量进行迭代，得到旋转增量及平移增量，并更新得到旋转分量和平移分量。一实施例中，为保证计算得到的数据的准确性，选取大量的基准点进行旋转增量及平移增量的计算，该基准点的选取数量优选为15个。将所选取的基准点通过牛顿法进行迭代计算出该次迭代的旋转增量和平移增量，以不断减小残差值，有效提高鲁棒性。在迭代过程中，基准点的选取数量是在不断地递增的。

在采用牛顿法求解时，需要求出残差方程关于变量的一阶导数，即求出残差方程的雅可比矩阵（Jacobian矩阵），例如，基准点的选取数量为15个，则该残差方程组为30个，因此Jacobian矩阵的阶数为30×6，该Jacobian矩阵如以下公式所示：

i＝1...30,j＝1...6,p_j表示体模中的一基准点。由此得到透视方程：

$\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]=\left[\begin{array}{c}-f{X}_i/s_x{Z}_i+o_x\\ -f{Y}_i/s_y{Z}_i+o_y\end{array}\right]=\left[\begin{array}{c}-f(r_{11}{X}_i^m+r_{12}{Y}_i^m+r_{13}{Z}_i^m+t_1)/s_x(r_{31}{X}_i^m+r_{32}{Y}_i^m+r_{33}{Z}_i^m+t_3)+o_x\\ -f(r_{21}{X}_i^m+r_{22}{Y}_i^m+r_{23}{Z}_i^m+t_2)/s_y(r_{31}{X}_i^m+r_{32}{Y}_i^m+r_{33}{Z}_i^m+t_3)+o_y\end{array}\right]$

其中，f为焦距，s_x和s_y表示体素的物理尺寸，o_x和o_y为投影图像的中心坐标，可在基准点标注时获得，t_i为平移分量，r_ij为旋转分量，

为体模中基准点的世界坐标，(X_i,Y_i,Z_i)为体模中基准点的物理坐标。

关于平移变量的Jacobian矩阵的为：

$\frac{\&PartialD;x_i}{\&PartialD;t_1}=\frac{-f}{s_x{Z}_i}$ $\frac{\&PartialD;x_i}{\&PartialD;t_2}=0$ $\frac{\&PartialD;x_i}{\&PartialD;t_3}=\frac{{\mathrm{fX}}_i}{s_x{Z}_i^2}$

$\frac{{\&PartialD;y}_i}{\&PartialD;t_1}=0$ $\frac{\&PartialD;y_i}{\&PartialD;t_2}=\frac{-f}{s_y{Z}_i}$ $\frac{\&PartialD;y_i}{\&PartialD;t_3}=\frac{f{Y}_i}{s_y{Z}_i^2}$

假设

为模型空间中的一点，与P_i相对应。

$\frac{\&PartialD;P_i^m}{\&PartialD;{\mathrm{\&Phi;}}_1}={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T\&times;P_i^m$

$\frac{\&PartialD;P_i^m}{\&PartialD;{\mathrm{\&Phi;}}_2}={\left[\begin{array}{ccc}0& 1& 0\end{array}\right]}^T\&times;P_i^m$

$\frac{\&PartialD;P_i^m}{\&PartialD;{\mathrm{\&Phi;}}_3}={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T\&times;P_i^m$

其中，Φ_j为旋转分量。

由于到P_i为旋转变换，因此有：

因此，关于旋转变量的Jacobian矩阵的为：

$\frac{\&PartialD;x_i}{\&PartialD;{\mathrm{\&Phi;}}_j}=-\frac{f(Z_i\frac{\&PartialD;X_i}{\&PartialD;{\mathrm{\&Phi;}}_j}-X_i\frac{\&PartialD;Z_i}{\&PartialD;{\mathrm{\&Phi;}}_j})}{s_x{Z}_i^2}$ $\frac{\&PartialD;y_i}{\&PartialD;{\mathrm{\&Phi;}}_j}=-\frac{f(Z_i\frac{\&PartialD;Y_i}{\&PartialD;{\mathrm{\&Phi;}}_j}-Y_i\frac{\&PartialD;Z_i}{\&PartialD;{\mathrm{\&Phi;}}_j})}{s_y{Z}_i^2}$

由 $\begin{array}{ccc}\frac{\&PartialD;X_i}{\&PartialD;{\mathrm{\&Phi;}}_j}& \frac{\&PartialD;Y_i}{\&PartialD;{\mathrm{\&Phi;}}_j}& \frac{\&PartialD;Z_i}{\&PartialD;{\mathrm{\&Phi;}}_j}\end{array}$ 的求解可以得到的求出，即

$\frac{\&PartialD;\left(\mathrm{\&delta;}{x}_i\right)}{\&PartialD;p_j}=\frac{\&PartialD;x_i}{\&PartialD;p_j}$ $\frac{\&PartialD;\left(\mathrm{\&delta;}{y}_i\right)}{\&PartialD;p_j}=\frac{\&PartialD;y_i}{\&PartialD;p_j}$

此时，由牛顿法可知，建立残差方程的一阶逼近为：

写成方程组的形式则是：JΔp＝δ

其中，J为Jacobian矩阵。

由旋转增量ΔΦ_j和平移增量Δt_j分别更新旋转分量和平移分量，详细公式如下所示：

$1,\left[\begin{array}{cc}{R}_{k+1}& T_{k+1}\\ 0& 1\end{array}\right]=\left[\begin{array}{cc}{R}_k& T_k\\ 0& 1\end{array}\right]\left[\begin{array}{cc}\mathrm{\&Delta;}{R}_k& \mathrm{\&Delta;}{T}_k\\ 0& 1\end{array}\right]$

2、T_k+1＝T_k+ΔT_k,Φ_k+1＝Φ_k+ΔΦ_k

3、T_k+1＝T_k+ΔT_k,R_k+1＝R_kΔR_k

在步骤S146中，判断相邻的两次选取基准点所得到的旋转分量之差及平移分量之差是否小于预设值或迭代次数是否大于预设次数，是，则终止迭代，反之，则返回步骤S144。通过牛顿法的不断迭代，来不断地降低残差值δ，在理想情况下残差值δ=0。为减少计算量，在投影图像的质量要求较低的情况下，其迭代次数至少1次。

在步骤S20中，根据三维运动矩阵，对非等中心投影图像转化，得到等中心投影图像。一实施方式中，根据三维运动矩阵，来获知C型臂姿态，从而结合射线锥束投影的几何原理，将非等中心投影图像转化为等中心投影图像。

具体地，如图7所示，步骤S20的过程具体是：

在步骤S202中，通过投影关系得到投影图像的投影坐标，根据旋转分量和平移分量计算得到圆心偏移量，并根据圆心偏移量和锥束射线的扫描角度，还原投影图像的物理坐标。一实施例中，图8示出了非等中心投影的几何关系图，射线源s绕Z轴旋转，β为C型臂的扫描角度，α为射线相对于旋转轴Z轴的角度，物理坐标v(x,y,z)为空间中的一点，为得到等中心投影图像中体素的坐标v′(p,q)，在射线源s进行扫描的过程中，可在投影平面上遍历所有体素的投影，以得到其坐标v′(p,q)，大小为row*col的投影图像，其计算量为o(row*col)。根据投影平面中体素的坐标v″(p′,q′)、圆心偏移量 $\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ \mathrm{\&Delta;z}\end{array}\right]$ 及几何关系，得到如下计算公式，以解出v(x,y,z)。

$\left\{\begin{array}{c}{p}^{\&prime;}=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*(y\cos \mathrm{\&beta;}-x\sin \mathrm{\&beta;}-\mathrm{\&Delta;}y\cos \mathrm{\&beta;}+\mathrm{\&Delta;}x\sin \mathrm{\&beta;})\\ q^{\&prime;}=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*(z-\mathrm{\&Delta;z})\\ f(x,y,z)=0\end{array}\right.$

其中，d_SM＝d_So+d_oM，f(x,y,z)为平面方程。d_oM＝xcosβ+ysinβ-Δxcosβ-Δysinβ，d_so为射线源到旋转中心的距离，d_SD为定值。

在步骤S204中，转换物理坐标得到等中心的投影图像。一实施例中，根据以下公式计算得到等中心的投影图像中每一点坐标v′(p,q)。

$p=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*(y\cos \left(\mathrm{\&beta;}\right)-x\sin \left(\mathrm{\&beta;}\right))$

$q=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*z$

对于上述还原投影图像的物理坐标和等中心的投影图像的求解过程，现结合图8和图9详细阐述该推导过程。在进行锥束图像重建时，要求所对应的投影图像为等中心的，即C型臂在扫描物体时，必须保证射线源到物体中心的距离不变。但是由于设计工艺等原因，在进行C型臂的扫描时，所得到的投影图像不能保证是等中心的，因此，需要根据射线锥束投影的几何原理，将非等中心投影图像转化为等中心投影图像。等中心投影图像的几何关系如图9所示，其中，s为射线源，(p,q)为C型臂扫描物体，并在平板探测器上得到投影图像的坐标，v(x,y,z)为物体的物理坐标，d_SO为射线源到旋转中心的距离，d_OD为旋转中心到平板探测器的距离。射线源沿Z轴等中心圆周扫描，并沿中心射线建立坐标系st，则v(x,y,z)在st坐标***下的坐标为：

$\left[\begin{array}{ccc}\cos \mathrm{\&beta;}& \sin \mathrm{\&beta;}& 0\\ -\sin \mathrm{\&beta;}& \cos \mathrm{\&beta;}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}s\\ t\\ z\end{array}\right]$

线段VM垂直于s轴，即M的坐标为(xcosβ+ysinβ，0)，VN为平行于t轴，且垂直于sz平面，射线源中的中心射线的投影点为D，V点的投影点为V′，则需要找到V′在投影坐标***pq中的投影坐标v′(p,q)。

d_om＝xcosβ+ysinβ，由三角形相似，可以得出，

$d_{\mathrm{VM}}=\sqrt{x^2+y^2+z^2-d_{\mathrm{om}}^2},$ d_sM＝d_so+d_oM

而d_SD为定值，则p＝d_V′D×sinα，q＝d_V′D×cosα

同样由ΔVNM和投影平面平行，得

$\sin \left(a\right)=\frac{d_{\mathrm{vN}}}{d_{\mathrm{vm}}}=\frac{y\cos \left(\mathrm{\&beta;}\right)-x\sin \left(\mathrm{\&beta;}\right)}{d_{\mathrm{vm}}}$

$\cos \left(\mathrm{\&alpha;}\right)=\frac{d_{\mathrm{NM}}}{d_{\mathrm{VM}}}=\frac{z}{d_{\mathrm{VM}}}$

综合上述公式，可以得出

$p=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*(y\cos \left(\mathrm{\&beta;}\right)-x\sin \left(\mathrm{\&beta;}\right))$

$q=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*z$

非等中心投影图像主要是由于C型臂的旋转中心发生偏移，不再对应物体中心，如图8所示，在非等中心投影图像中，射线源仍绕Z轴旋转，β为扫描角度，α为射线相对于旋转轴Z轴的角度，s为射线源，圆心偏移量为 $\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ \mathrm{\&Delta;z}\end{array}\right],$ 物理坐标v(x,y,z)为空间中的一点，v(x,y,z)的投影点为v″(p′,q′)，同样如前所述的所推导的三角关系可以得到：

p′＝d_V′D×sinαq′＝d_V′D×cosα

$\sin \mathrm{\&alpha;}=\frac{d_{\mathrm{VN}}}{d_{\mathrm{VM}}}$ $\cos \mathrm{\&alpha;}=\frac{d_{\mathrm{NM}}}{d_{\mathrm{VM}}}$

由几何关系得到M点在st坐标下的位置，即M(xcosβ+ysinβ-Δxcosβ-Δysinβ，0)

d_vn＝xcosβ+ysinβ-(Δxcosβ+Δysinβ)

d_MN＝z-Δz

$d_{\mathrm{VM}}=\sqrt{{(x-\mathrm{\&Delta;x})}^2+{(y-\mathrm{\&Delta;y})}^2+{(z-\mathrm{\&Delta;z})}^2-d_{\mathrm{OM}}^2}$

综上：

$p^{\&prime;}=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}\&times;(y\cos \mathrm{\&beta;}-x\sin \mathrm{\&beta;}-(\mathrm{\&Delta;}y\cos \mathrm{\&beta;}-\mathrm{\&Delta;}x\sin \mathrm{\&beta;}))$

$q^{\&prime;}=\frac{d_{\mathrm{SM}}}{d_{\mathrm{SM}}}*(z-\mathrm{\&Delta;z})$

d_SM＝d_So+d_oM，d_So为一定值。

d_OM＝xcosβ+ysinβ-Δxcosβ+Δysinβ

当v(x,y,z)、旋转角度β和圆心偏移量 $\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ \mathrm{\&Delta;z}\end{array}\right]$ 已知，即可分别求出其在等中心图像和非等中心图像下的坐标，反之，若已知其在非等中心图像下的坐标v″(p′,q′)、旋转角度β和圆心偏移量 $\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ \mathrm{\&Delta;z}\end{array}\right],$ 则找到一组v(x,y,z)，该v(x,y,z)满足非等中心图像下的坐标为v″(p′,q′)，就可求得v′(p,q)。

在步骤S30中，重建锥束图像。一实施方式中，为适应有限角投影、噪声状况差的情况，将投影图像重建的过程归结为计算每个体素的衰减系数的过程。具体地，如图10所示，重建锥束图像的过程具体是：

在步骤S302中，通过等中心投影图像中体素的衰减系数和投影图像获得体素值向量。一实施例中，通过公式W^θv＝P^θ计算得到体素值向量v＝[v₁ v₂...v_N]^T，体素值向量包括了N＝N₁×N₂×N₃个元素，这些体素构成了N₁×N₂×N₃像素大小的目标重建体，p^θ＝[p₁...p_m*m]^T是投影数据向量，若共有k个投影角度，则每一个投影角度的投影数据由L（m×n）个像素组成，则R=m×n。投影数据值p_j ^θ为像素j在投影角度为θ时获得，则在θ角度有m×m条射线，W^θ是一个R×N大小的权重矩阵，记W_i，j ^θ为在θ投影角度下从体素v_i到像素j的衰减系数，在不同角度下所对应的值是不同的，i＝1...N，j＝1...R。由离散锥束投影值，在θ角度下的投影射线为有限条，可以对同一角度下的不同射线采用并行计算加速，即采用显卡GPU加速，迭代重建采用SART算法。

在步骤S304中，修正等中心投影图像，得到修正值。一实施例中，通过以下公式得到等中心投影图像的修正值。

${{\mathrm{\&Delta;p}}_j}^{\mathrm{\&theta;}}=p_j^{\mathrm{\&theta;}}-{p^{*\mathrm{\&theta;}}}_j$

其中，为θ角度时j射线的等中心投影图像值。

在步骤S306中，根据修正值和体素的衰减系数进行反投影，更新体素值向量。一实施例中，通过公式

进行反投影，并通过公式

更新体素v_i值。其中，α_jθ为θ角度下穿过第j条射线的所有体素的衰减系数的平方和，

γ为松弛算子，一般选取（0~2），γ是在上述投影反投影的过程中，进行逐角度更新体素，直到判断体素收敛而得到的，即收敛次数大于预设次数或者更新后的体素投影与真投影之间的误差小于预设值。

在进行锥束图像的重建过程中，最为耗时的投影和反投影过程中，每个过程都可分为基于像素或者体素的单独计算过程，可用基于图形处理器（GraphicProcessing Unit，简称GPU）并行计算的方法加速。

图11示出了一实施例中的锥束图像重建装置，包括跟踪模块10、中心变换模块20以及重建模块30，其中：

跟踪模块10用于扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵。一实施方式中，跟踪模块10对C型臂的运动姿态进行跟踪，以修正非中心投影，具体地，跟踪模块10至少包括：

提取单元102，用于提取体模基准点。如前所述，体模为镶嵌预设数量铅球的物体，该铅球呈螺旋形排列，因此在其投影上表现为若干像素的阴影。体模中的铅球质地均匀，内部没有瑕疵，半径在1.5~3mm内，间距为3~8mm，从而保证铅球的投影不会产生重叠。一实施例中，提取单元102对投影图像进行高斯滤波后进行二值化处理，并通过连通成分标记进行基准点的标注。该提取单元102通过图像二值化的方法来进行基准点的分割，假定原灰度图像为f(x,y)，则二值化后的新图像为

$p(x,y)=\left\{\begin{array}{c}1,f(x,y)\&GreaterEqual;T\\ 0,f(x,y)<T\end{array}\right.$

其中，T为f(x,y)领域的平均值，

其中，k为阈值，M为k邻域内像素个数。

图像二值化之后，提取单元102进行逐个体素扫描图像，扫描方式是沿图像列顺序遍历，然后沿图像行顺序遍历。查找当前体素的相邻体素，如果体素不是背景且没有相邻体素，则标识该体素，并继续遍历，如果体素不是背景并且有相邻的体素，该体素和相邻的体素都分配了不同的标识数，则将该体素与其相邻的体素都标示为同一连通成分，即把该体素与其相邻体素的所有标记数中的最小值分配给这一连通成分，并存储其相邻关系及标识数，最后用最小的标识数重新标识体素。通过连通成分标记既可以得到所有的连通成分，又可以在此基础上独立计算连通区域的统计特性。由于体模中的基准点在投影图像上为若干像素的集合，因此对每一基准点均用像素的形心(i,j)来表示，具体公式如下所示：

$\stackrel{\&OverBar;}{i}=\frac{\underset{i=r}{\overset{m-1}{\mathrm{\&Sigma;}}}i\&times;f(x,y)}{\underset{i=r}{\overset{m-1}{\mathrm{\&Sigma;}}}f(x,y)}$ $\stackrel{\&OverBar;}{j}=\frac{\underset{j=c}{\overset{n-1}{\mathrm{\&Sigma;}}}j\&times;f(x,y)}{\underset{j=c}{\overset{n-1}{\mathrm{\&Sigma;}}}f(x,y)}$

其中i，j所选用的区域为基准点所在的联通区域。

在其他实施例中，提取单元102剔除投影图像中CT值大于第一阈值的区域，得到初选的投影区域，在该初选的投影区域中提取CT值大于第二阈值的感兴趣的投影区域，第一阈值大于第二阈值，并通过连通成分标记进行基准点的标注。第一阈值优选为3000Hounsfield uints，第二阈值优选为2000Hounsfield uints。

估计单元104，用于根据体模基准点，进行C型臂的运动姿态估计，得到三维运动矩阵。一实施例中，估计单元104选取基准点，设定三维运动矩阵初值，建立残差方程计算得到基准点初始旋转分量和初始平移分量，然后再次选取基准点，将初始的旋转分量和平移分量进行迭代，得到旋转增量及平移增量，更新旋转分量和平移分量，判断相邻的两次选取基准点所得到的旋转分量之差及平移分量之差是否小于预设阈值或迭代次数是否大于预设次数，是则终止迭代，反之，则再次选取基准点进行迭代，最后通过旋转增量及平移增量分别，其迭代次数至少1次。

在优选的实施例中，体模中设定了15个基准点，每间隔两个基准点选取一个，共选取5个基准点。三维运动矩阵为 $F_x=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& r_{14}\\ r_{21}& r_{22}& r_{23}& r_{24}\\ r_{31}& r_{32}& r_{33}& r_{34}\\ r_{41}& r_{42}& r_{43}& r_{44}\end{array}\right]$

其中，r_ij为旋转分量，T_ij为平移分量，F_x为变换矩阵，表示物体相对于射线源的变换，三维运动矩阵的假定初值为F_x＝[1000 0100 0010 0001]。

建立残差方程，并将基准点的投影坐标及标注导入残差方程中，得到初始旋转分量和初始平移分量，残差方程如下所示：

$\mathrm{\&delta;}{x}_i=x_i(R_k,T_k)-{\stackrel{\&OverBar;}{x}}_i$

$\mathrm{\&delta;}{y}_i=y_i(R_k,T_k)-{\stackrel{\&OverBar;}{y}}_i$

其中，(x_i,y_i)为基准点的投影坐标，

为基准点的标注，(δx_i,δy_i)为残差值，R_k为F_x的旋转矩阵，T_k为F_x的平移向量。

再次选取15个基准点，则该残差方程组为30个，因此Jacobian矩阵的阶数为30×6，该Jacobian矩阵如以下公式所示：

i＝1...30,j＝1...6,p_j表示体模中的一基准点。

由此得到透视方程：

$\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]=\left[\begin{array}{c}-f{X}_i/s_x{Z}_i+o_x\\ -f{Y}_i/s_y{Z}_i+o_y\end{array}\right]=\left[\begin{array}{c}-f(r_{11}{X}_i^m+r_{12}{Y}_i^m+r_{13}{Z}_i^m+t_1)/s_x(r_{31}{X}_i^m+r_{32}{Y}_i^m+r_{33}{Z}_i^m+t_3)+o_x\\ -f(r_{21}{X}_i^m+r_{22}{Y}_i^m+r_{23}{Z}_i^m+t_2)/s_y(r_{31}{X}_i^m+r_{32}{Y}_i^m+r_{33}{Z}_i^m+t_3)+o_y\end{array}\right]$

其中，f为焦距，s_x和s_y表示体素的物理尺寸，o_x和o_y为投影图像的中心坐标，可在基准点标注时获得，t_i为平移分量，r_ij为旋转分量，

为体模中基准点的世界坐标，(X_i,Y_i,Z_i)为体模中基准点的物理坐标。

关于平移变量的Jacobian矩阵的为：

$\frac{\&PartialD;x_i}{\&PartialD;t_1}=\frac{-f}{s_x{Z}_i}$ $\frac{\&PartialD;x_i}{\&PartialD;t_2}=0$ $\frac{\&PartialD;x_i}{\&PartialD;t_3}=\frac{{\mathrm{fX}}_i}{s_x{Z}_i^2}$

$\frac{{\&PartialD;y}_i}{\&PartialD;t_1}=0$ $\frac{\&PartialD;y_i}{\&PartialD;t_2}=\frac{-f}{s_y{Z}_i}$ $\frac{\&PartialD;y_i}{\&PartialD;t_3}=\frac{f{Y}_i}{s_y{Z}_i^2}$

假设

为模型空间中的一点，与P_i相对应。

$\frac{\&PartialD;P_i^m}{\&PartialD;{\mathrm{\&Phi;}}_1}={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T\&times;P_i^m$

$\frac{\&PartialD;P_i^m}{\&PartialD;{\mathrm{\&Phi;}}_2}={\left[\begin{array}{ccc}0& 1& 0\end{array}\right]}^T\&times;P_i^m$

$\frac{\&PartialD;P_i^m}{\&PartialD;{\mathrm{\&Phi;}}_3}={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T\&times;P_i^m$

其中，Φ_j为旋转分量。

由于到P_i为旋转变换，因此有：

因此，关于旋转变量的Jacobian矩阵的为：

$\frac{\&PartialD;x_i}{\&PartialD;{\mathrm{\&Phi;}}_j}=-\frac{f(Z_i\frac{\&PartialD;X_i}{\&PartialD;{\mathrm{\&Phi;}}_j}-X_i\frac{\&PartialD;Z_i}{\&PartialD;{\mathrm{\&Phi;}}_j})}{s_x{Z}_i^2}$ $\frac{\&PartialD;y_i}{\&PartialD;{\mathrm{\&Phi;}}_j}=-\frac{f(Z_i\frac{\&PartialD;Y_i}{\&PartialD;{\mathrm{\&Phi;}}_j}-Y_i\frac{\&PartialD;Z_i}{\&PartialD;{\mathrm{\&Phi;}}_j})}{s_y{Z}_i^2}$

由 $\begin{array}{ccc}\frac{\&PartialD;X_i}{\&PartialD;{\mathrm{\&Phi;}}_j}& \frac{\&PartialD;Y_i}{\&PartialD;{\mathrm{\&Phi;}}_j}& \frac{\&PartialD;Z_i}{\&PartialD;{\mathrm{\&Phi;}}_j}\end{array}$ 的求解可以得到

的求出，即

$\frac{\&PartialD;\left(\mathrm{\&delta;}{x}_i\right)}{\&PartialD;p_j}=\frac{\&PartialD;x_i}{\&PartialD;p_j}$ $\frac{\&PartialD;\left(\mathrm{\&delta;}{y}_i\right)}{\&PartialD;p_j}=\frac{\&PartialD;y_i}{\&PartialD;p_j}$

此时，由牛顿法可知，建立残差方程的一阶逼近为：

写成方程组的形式则是：JΔp＝δ

其中，J为Jacobian矩阵。

由旋转增量ΔΦ_j和平移增量Δt_j分别更新旋转分量和平移分量，详细公式如下所示：

$1,\left[\begin{array}{cc}{R}_{k+1}& T_{k+1}\\ 0& 1\end{array}\right]=\left[\begin{array}{cc}{R}_k& T_k\\ 0& 1\end{array}\right]\left[\begin{array}{cc}\mathrm{\&Delta;}{R}_k& \mathrm{\&Delta;}{T}_k\\ 0& 1\end{array}\right]$

2、T_k+1＝T_k+ΔT_k,Φ_k+1＝Φ_k+ΔΦ_k

3、T_k+1＝T_k+ΔT_k,R_k+1＝R_kΔR_k

中心变换模块20，用于根据三维运动矩阵，对非等中心投影图像转化，得到等中心投影的投影图像。具体地，中心变换模块20包括：

还原单元202，用于通过投影关系得到非等中心投影图像的投影坐标，根据旋转分量和平移分量计算得到圆心偏移量，并根据圆心偏移量和锥束射线的扫描角度，还原投影图像的物理坐标。一实施例中，还原单元202通过如下计算公式，以解出投影图像的物理坐标v(x,y,z)。

$\left\{\begin{array}{c}{p}^{\&prime;}=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*(y\cos \mathrm{\&beta;}-x\sin \mathrm{\&beta;}-\mathrm{\&Delta;}y\cos \mathrm{\&beta;}+\mathrm{\&Delta;}x\sin \mathrm{\&beta;})\\ q^{\&prime;}=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*(z-\mathrm{\&Delta;z})\\ f(x,y,z)=0\end{array}\right.$

其中，β为C型臂的扫描角度， $\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ \mathrm{\&Delta;z}\end{array}\right]$ 为圆心偏移量，d_SM＝d_So+d_oM，f(x,y,z)为平面方程，v″(p′,q′)为投影平面中体素的坐标。d_oM＝xcosβ+ysinβ-Δxcosβ-Δysinβ，d_so为射线源到旋转中心的距离，d_SD为定值。

转换单元204，用于转换物理坐标得到等中心的投影图像。该转换单元204根据以下公式计算得到等中心的投影图像中每一点坐标v′(p,q)。

$p=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*(y\cos \left(\mathrm{\&beta;}\right)-x\sin \left(\mathrm{\&beta;}\right))$

$q=\frac{d_{\mathrm{SD}}}{d_{\mathrm{SM}}}*z$

重建模块30，用于重建锥束图像。一实施方式中，重建模块30将投影图像重建的过程归结为计算每个体素的衰减系数的过程。具体地，重建模块30包括：

投影单元302，用于通过等中心投影图像中体素的衰减系数和投影图像获得体素值向量。一实施例中，投影单元302通过公式W^θv＝P^θ计算得到体素值向量v＝[v₁ v₂...v_N]^T，体素值向量包括了N＝N₁×N₂×N₃个元素，这些体素构成了N₁×N₂×N₃体素大小的目标重建体，p^θ＝[p₁...p_m*m]^T是投影数据向量，若共有k个投影角度，则每一个投影角度的投影数据由L（m×n）个像素组成，则R=m×n。投影数据值p_j ^θ都是在投影角度为θ时获得，则在θ角度有m×n条射线，W^θ是一个R×N大小的权重矩阵，记W_i,j ^θ为在投影角度θ从体素v_i到像素j的衰减系数，在不同角度下所对应的值是不同的，i＝1...N，j＝1...R。

修正单元304，用于修正等中心投影图像，得到修正值。一实施例中，修正单元304通过以下公式得到等中心投影图像的修正值，图12示出了一实施例的原始等中心投影图像与等中心图像之差。

${{\mathrm{\&Delta;p}}_j}^{\mathrm{\&theta;}}=p_j^{\mathrm{\&theta;}}-{p^{*\mathrm{\&theta;}}}_j$

其中，

为θ角度时获得的原始投影图像。

反投影单元306，用于根据修正值和体素的衰减系数进行反投影，更新体素值向量。一实施例中，反投影单元306通过公式

进行反投影，并通过公式

更新体素值向量。其中，α_jθ为θ角度下穿过第j条射线的所有体素的衰减系数的平方和

γ为松弛算子，一般选取（0~2），γ是在上述投影反投影的过程中，进行逐角度更新体素，直到判断体素收敛而得到的，即收敛次数大于预设次数或者更新后的体素投影与真投影之间的误差小于预设值。

请结合参阅图13至图16c为通过上述锥束图像重建方法及装置得到的实例。如图13与14所示，图13比图14更为平滑，成像质量更高。分别将图15a及图16a、图15b及图16b、图15c及图16c进行比较，可以清楚地看出图16a、图16b和图16c这一组脊骨三维截面图与图15a、图15b和图15c的这一组脊骨原始三维截面图非常接近，由上述锥束图像重建方法及装置得到的脊骨原始三维截面图成像质量非常高。

上述锥束图像重建方法及装置通过跟踪C型臂的运动姿态，将扫描得到的非等中心投影图像结合投影的几何关系，在不增加额外跟踪设备的基础上，转化为等中心图像，有效地降低了成本，提高投影图像的准确性，提高了投影图像的成像质量及其鲁棒性。

上述锥束图像重建方法及装置通过多次运动姿态估计，来得到三维运动矩阵，收敛性好，有效地降低了投影点与基准点坐标之间的误差。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (12)

1.一种锥束图像重建方法，包括如下步骤：

扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵；

所述扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵的具体过程为：

提取体模基准点，根据所述体模基准点进行C型臂的运动姿态估计，运用投影图像上至少3组基准点建立方程组得到三维运动矩阵；

所述姿态估计的步骤是：

选取基准点，设定三维运动矩阵初值，建立残差方程计算得到所述基准点的初始旋转分量和初始平移分量，残差方程如下所示：

$\mathrm{\&delta;}{x}_i=x_i(R_k,T_k)-{\stackrel{\&OverBar;}{x}}_i$

$\mathrm{\&delta;}{y}_i=y_i(R_k,T_k)-{\stackrel{\&OverBar;}{y}}_i$

其中，(x_i,y_i)为基准点的投影坐标，

为基准点的标注，(δx_i,δy_i)为残差值，R_k为F_x的旋转矩阵，T_k为F_x的平移向量，F_x为三维运动矩阵；

再次选取基准点，将初始的旋转分量和平移分量进行迭代，得到旋转增量及平移增量，并更新得到旋转分量和平移分量；

判断相邻的两次选取基准点所得到的旋转分量之差及平移分量之差是否小于预设阈值或迭代次数是否大于预设次数，是，则终止迭代，反之，则再次选取基准点进行迭代；

根据所述三维运动矩阵，对非等中心投影图像转化，得到等中心投影图像；

根据所述等中心投影图像，重建锥束图像。

2.根据权利要求1所述的锥束图像重建方法，其特征在于，所述提取体模基准点的步骤是：

对所述投影图像进行高斯滤波；

将滤波后得到的投影图像进行二值化处理；

通过连通成分标记进行基准点的标注。

3.根据权利要求1所述的锥束图像重建方法，其特征在于，所述提取体模基准点的步骤是：

剔除投影图像中CT值大于第一阈值的区域，得到初选的投影区域；

在所述初选的投影区域中提取CT值大于第二阈值的感兴趣的投影区域，所述第一阈值大于第二阈值；

通过连通成分标记进行基准点的标注。

4.根据权利要求1所述的锥束图像重建方法，其特征在于，所述迭代次数至少1次。

5.根据权利要求1所述的锥束图像重建方法，其特征在于，所述根据所述三维运动矩阵，对非等中心投影图像转化，得到等中心投影图像的步骤是：

通过投影关系得到投影图像的投影坐标，根据旋转分量和平移分量计算得到圆心偏移量，并根据所述圆心偏移量和锥束射线的扫描角度，还原所述投影图像的物理坐标；

转换所述物理坐标得到等中心的投影图像。

6.根据权利要求1所述的锥束图像重建方法，其特征在于，所述根据所述等中心投影图像，重建锥束图像的步骤是：

通过等中心投影图像中体素的衰减系数和所述投影图像获得体素值向量；

修正所述等中心投影图像，得到修正值；

根据所述修正值和体素的衰减系数进行反投影，更新体素值向量。

7.一种锥束图像重建装置，其特征在于，至少包括：

跟踪模块，用于扫描得到投影图像，并跟踪C型臂的运动姿态，计算得到三维运动矩阵；

所述跟踪模块包括提取单元和估计单元，所述提取单元用于提取体模基准点，所述估计单元用于根据所述体模基准点进行C型臂的运动姿态估计，运用投影图像上至少3组基准点建立方程组得到三维运动矩阵；

进一步的，所述跟踪模块选取基准点，设定三维运动矩阵初值，建立残差方程计算得到所述基准点的初始旋转分量和初始平移分量，残差方程如下所示：

$\mathrm{\&delta;}{x}_i=x_i(R_k,T_k)-{\stackrel{\&OverBar;}{x}}_i$

$\mathrm{\&delta;}{y}_i=y_i(R_k,T_k)-{\stackrel{\&OverBar;}{y}}_i$

其中，(x_i,y_i)为基准点的投影坐标，为基准点的标注，(δx_i,δy_i)为残差值，R_k为F_x的旋转矩阵，T_k为F_x的平移向量，F_x为三维运动矩阵；再次选取基准点，将初始的旋转分量和平移分量进行迭代，得到旋转增量及平移增量，并更新得到旋转分量和平移分量；判断相邻的两次选取基准点所得到的旋转分量之差及平移分量之差是否小于预设阈值或迭代次数是否大于预设次数，是，则终止迭代，反之，则再次选取基准点进行迭代；

中心变换模块，用于根据所述三维运动矩阵，对非等中心投影图像转化，得到等中心投影的投影图像；

重建模块，用于根据所述等中心投影图像，重建锥束图像。

8.根据权利要求7所述的锥束图像重建装置，其特征在于，所述提取单元对所述投影图像在高斯滤波后进行二值化处理，并通过连通成分标记进行基准点的标注。

9.根据权利要求7所述的锥束图像重建装置，其特征在于，所述提取单元剔除投影图像中CT值大于第一阈值的区域，得到初选的投影区域，在所述初选的投影区域中提取CT值大于第二阈值的感兴趣的投影区域，所述第一阈值大于第二阈值，并通过连通成分标记进行基准点的标注。

10.根据权利要求7所述的锥束图像重建装置，其特征在于，所述估计单元中的迭代次数至少1次。

11.根据权利要求7所述的锥束图像重建装置，其特征在于，所述中心变换模块包括：

还原单元，用于通过投影关系得到非等中心投影图像的投影坐标，根据旋转分量和平移分量计算得到圆心偏移量，并根据圆心偏移量和锥束射线的扫描角度，还原投影图像的物理坐标；

转换单元，用于转换所述物理坐标得到等中心的投影图像。

12.根据权利要求7所述的锥束图像重建装置，其特征在于，所述重建模块包括：

投影单元，用于通过等中心投影图像中体素的衰减系数和所述投影图像获得体素值向量；

修正单元，用于修正所述等中心投影图像，得到修正值；

反投影单元，用于根据所述修正值和体素的衰减系数进行反投影，更新体素值向量值。

Images