CN103258335A - 基于普通c臂机的脊柱横断面图像合成方法 - Google Patents

Abstract

一种基于普通C臂机的脊柱横断面图像合成方法，包括以下步骤：步骤一、对普通C臂机获得的2D锥束投影图像进行预处理，得到1D扇束投影数据；步骤二、对步骤一得到的扇束进行虚拟旋转再处理，计算得到虚拟扇束投影数据；步骤三、对2D断层图像合成方法中的相关参数进行调整，采用调整后的方法基于虚拟扇束投影数据合成脊柱横断面图像。本发明合成方法，原理简单，执行方便，除C臂机本身外无需购置或定制其他额外装置，可用于合成在临床诊断、手术规划、引导手术操作方面都有重要意义的脊柱横断面图像。

Description

基于普通C臂机的脊柱横断面图像合成方法

技术领域

本发明涉及计算机成像技术，特别涉及一种基于普通C臂机的脊柱横断面图像合成方法。

背景技术

在医院外科临床应用中，尤其是脊柱外科领域，C臂机是一种常用的设备。借助C臂机获取的图像，医生可对患者脊柱的病损状况进行判断，还可对患者进行手术规划和植入锥弓根螺钉等手术操作。但是，通过普通C臂机仅可获取患者的投影图像，而临床上对病情的判断及引导手术操作意义重大的脊柱横断面图像却不能通过C臂机直接获取。

计算机断层成像技术利用X线对患者进行扫描获取投影图像序列，利用特定的算法以投影图像序列为基础合成人体横断面图像。经对现有技术的文献检索发现，Feldkamp等人提出的FDK算法(Practical cone-beam reconstruction，J.Opt.Soc.Am.A1，612-619，1984)适用基于扇束和锥束投影的断层图像合成，但其要求具备360°范围内全扫描路径投影数据，而普通C臂旋转扫描角度范围一般不超过180°，不能直接采用FDK算法。Noo和Kudo等人提出的超短扫描路径断层图像合成算法(Image reconstruction from fan-beam projections on less than ashort scan，J.Phys.Med.Biol.47，2525-2546，2002；New Super-Short-ScanAlgorithms for Fan-Beam and Cone-Beam Reconstruction，IEEE Nuclear Science andMedical Imaging Symposium，Record，Norfolk，USA，2003，2：902-906)允许扫描范围低于180°，但其要求设备必须满足等中心旋转扫描条件，而普通C臂锥束中轴线一般不通过C臂旋转中心，其并不满足等中心旋转扫描条件，也不能直接采用这些算法。

发明内容

本发明的目的，就在于克服上述已有技术的不足，提供一种基于普通C臂机的脊柱横断面图像合成方法，以便依托目前国内大多数医院都具备的普通C臂机合成在临床上意义重大的横断面图像。

为实现上述目的，本发明采用了以下技术方案：一种基于普通C臂机的脊柱横断面图像合成方法，包括以下步骤：

步骤一、对普通C臂机获得的2D锥束投影图像进行预处理，得到1D扇束投影数据；

步骤二、对步骤一得到的扇束进行虚拟旋转再处理，计算得到虚拟扇束投影数据；

步骤三、对2D断层图像合成方法中的相关参数进行调整，采用调整后的方法基于虚拟扇束投影数据合成脊柱横断面图像。

步骤一中所述的预处理，是指首先对2D锥束投影图像进行中值滤波除噪，然后提取出图像中间的3列或2列像素的灰度值，对提取的像素的灰度值取均值构成1D扇束投影数据；

所述1D扇束投影与普通C臂机获得的2D锥束投影的关系为：扇束光源点与C臂X线锥束光源点重合，扇束中轴线与C臂X线锥束中轴线重合，扇束光源点到其成像直线之间的距离等于C臂X线锥束光源点到其成像平面之间的距离，扇束成像直线位于C臂X线锥束成像平面上，扇束的扇形角等于C臂X线锥束的锥角γ_m。

步骤二中所述的再处理，是指对步骤一所得扇束绕其光源点进行虚拟旋转得到虚拟扇束，通过计算后得到对应的虚拟扇束投影数据；

所述虚拟扇束与扇束的关系为：虚拟扇束光源点到其成像直线之间的距离等于扇束光源点到其成像直线之间的距离D，虚拟扇束光源点到其成像直线的垂线通过C臂旋转中心；对虚拟扇束投影和扇束投影分别建立笛卡尔坐标系，令两个坐标系的原点分别位于光源点在各自成像直线上的垂足，且令虚拟扇束投影坐标轴正向与扇束投影坐标轴正向所夹角度为θ，则两个坐标系下投影像素坐标及灰度值对应关系为：

$u=D.\tan (\arctan \left(\frac{u_v}{D}\right)-\mathrm{\θ}),$ $g_v(\mathrm{\λ},u_v)=g_c(\mathrm{\λ},u).\left(\frac{\sqrt{D^2+u_v^2}}{\sqrt{D^2+u^2}}\right)$

上式中，u为扇束投影在其笛卡尔坐标系中的像素坐标，u_v为虚拟扇束投影在其笛卡尔坐标系中的像素坐标，g_c(λ，u)为扇束投影在对应坐标处的像素灰度值，g_v(λ，u_v)为虚拟扇束投影在对应坐标处的像素灰度值，λ为C臂旋转角度。

步骤三中所述的对2D断层图像合成方法中的相关参数进行调整，具体包括：将光源点沿圆弧路径运动的半径R调整为光源点到C臂旋转中心的距离，将对扇束投影数据进行滤波处理的区间下限u_l调整为-D.tan(γ_m-θ)，上限u_u调整为D.tan(γ_m+θ)，上述各式中，γ_m为扇束的扇形角，D为扇束光源点到其成像直线之间的距离，θ为虚拟扇束投影坐标轴正向与扇束投影坐标轴正向所夹角度。

所述提取出图像中间的3列或2列像素的灰度值，是指当图像像素列数n为奇数时，提取出图像中间的3列像素的灰度值，当图像像素列数n为偶数时，提取出图像中间的2列像素的灰度值。

本发明的有益效果是，采用该方法，依托目前国内大多数医院都具备的普通非等中心旋转扫描C臂机，可合成在临床诊断、手术规划、引导手术操作方面都有重要意义的脊柱横断面图像。

本发明基于普通C臂机的脊柱横断面图像合成方法原理简单，执行方便，除C臂机本身外无需购置或定制其他额外装置。

附图说明

图1为本发明基于普通C臂机的脊柱横断面图像合成方法的流程图；

图2为本发明步骤一中锥束投影与扇束投影的关系示意图；

图3为本发明步骤二中扇束与虚拟扇束的关系示意图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

参见图1，本发明基于普通C臂机的脊柱横断面图像合成方法包含三个步骤：

第1步，对普通C臂2D锥束投影图像进行预处理，得到1D扇束投影数据。首先对C臂在每个旋转角度λ下获取的2D锥束投影图像g_m(λ，u，v)进行中值滤波除噪(u，v分别为图像中像素的行、列坐标)，得到除噪后投影图像g_e(λ，u，v)。然后再从图像g_e(λ，u，v)中提取出其中间3列像素的灰度值(若图像像素列数n为奇数)：g_c1(λ，u)＝g_e(λ，u，v_(((n+1)/2)-1)，g_c2(λ，u)＝g_e(λ，u，v_((n+1)/2))，g_c3(λ，u)＝g_e(λ，u，v_(((n+1)/2)+1)，或中间2列像素的灰度值(若图像像素列数n为偶数)：g_cl(λ，u)＝g_e(λ，u，v_(n/2))，g_c2(λ，u)＝g_e(λ，u，v_((n/2)+1))。最后对中间3列(或2列)像素的灰度值取均值构成1D扇束投影数据g_c(λ，u)。配合参见图2，所得扇束光源点21与C臂X线锥束1的光源点重合，所得扇束中轴线23与C臂X线锥束中轴线重合，所得扇束光源点21到其成像直线25之间的距离等于C臂X线锥束光源点到其成像平面15之间的距离，所得扇束成像直线25位于C臂X线锥束成像平面15上，所得扇束的扇形角22等于C臂X线锥束的锥角，设其为γ_m。

第2步，配合参见图2和图3，对扇束2进行虚拟旋转再处理，得到虚拟扇束3，并计算得到虚拟扇束投影数据。虚拟扇束光源点与扇束光源点21重合，虚拟扇束光源点到其成像直线35之间的距离等于扇束光源点21到其成像直线25之间的距离，设其为D，虚拟扇束光源点到其成像直线35的垂线33通过C臂旋转中心4。设θ为虚拟扇束投影成像直线35之垂线与扇束投影成像直线25之垂线所夹锐角。对虚拟扇束投影建立笛卡尔坐标系，令原点位于其光源点在其成像直线35上的垂足34。对扇束投影建立笛卡尔坐标系，令原点位于其光源点21在其成像直线25上的垂足24。令虚拟扇束3之投影坐标轴正向与扇束2之投影坐标轴正向所夹角度为θ，令两个坐标系下投影像素坐标及灰度值对应关系为：

$u=D.\tan (\arctan \left(\frac{u_v}{D}\right)-\mathrm{\θ}),$ $g_v(\mathrm{\λ},u_v)=g_c(\mathrm{\λ},u).\left(\frac{\sqrt{D^2+u_v^2}}{\sqrt{D^2+u^2}}\right)$

上式中，u和u_v分别为扇束投影及虚拟扇束投影在各自坐标系中的像素坐标，g(u)和g_v(u_v)分别为扇束投影及虚拟扇束投影在对应坐标处的像素灰度值，λ为C臂旋转角度。

第3步，对2D断层图像合成方法中的相关参数进行调整，采用调整后的方法基于虚拟扇束投影数据合成横断面图像。配合参见图2和图3，将光源点沿圆弧路径运动的半径R调整为光源点21到C臂旋转中心4之间的距离，将扇束投影数据滤波处理区间的下限u_l调整为-D.tan(γ_m-θ)，上限u_u调整为D.tan(γ_m+θ)。采用公式 $f\left(\stackrel{\→}{x}\right)=\frac{1}{2\mathrm{\π}}{\∫}_0^{{\mathrm{\λ}}_u}\mathrm{d\λ}\frac{1}{R+\stackrel{\→}{x}.{\stackrel{\→}{e}}_1}{\left[w(\mathrm{\λ},\tilde{u})g_f(\mathrm{\λ},\tilde{u})\right]}_{\tilde{u}=\tilde{u}*(\mathrm{\λ},\stackrel{\→}{x})},$ 基于虚拟扇束投影数据合成横断面图像。和

分别表示横断面图像像素坐标及像素灰度值。γ_m的含义见步骤1，D和θ的含义见步骤2。其他中间计算公式具体如下：

$g_f(\mathrm{\λ},\tilde{u})={\∫}_{u_l}^{u_u}{\mathrm{duh}}_H(\tilde{u}-u_v)\frac{D}{\sqrt{D^2+{u_v}^2}}(\frac{\∂}{\∂\mathrm{\λ}}+\frac{D^2+{u_v}^2}{D}\frac{\∂}{\∂u_v})g_v(\mathrm{\λ},u_v),$

$u_l=-D.\tan ({\mathrm{\γ}}_m-\mathrm{\θ}),u_u=D.\tan ({\mathrm{\γ}}_m+\mathrm{\θ}),h_H(\tilde{u}-u)=\frac{1}{\mathrm{\π}(\tilde{u}-u)},$

$\frac{\∂}{\∂\mathrm{\λ}}{g}_v(\mathrm{\λ},u_v)\≈\frac{g_v(\mathrm{\λ}+\mathrm{\Δ\λ},u_v)-g_v(\mathrm{\λ}-\mathrm{\Δ\λ},u_v)}{2\mathrm{\Δ\λ}},$

$\frac{\∂}{\∂u_v}{g}_v(\mathrm{\λ},u_v)\≈\frac{g_v(\mathrm{\λ},u_v+\mathrm{\Δ}{u}_v)-g_v(\mathrm{\λ},u_v-\mathrm{\Δ}{u}_v)}{2\mathrm{\Δ}{u}_v},$

$f\left(\stackrel{\→}{x}\right)=\frac{1}{2\mathrm{\π}}{\∫}_0^{{\mathrm{\λ}}_u}\mathrm{d\λ}\frac{1}{R+\stackrel{\→}{x}.{\stackrel{\→}{e}}_1}{\left[w(\mathrm{\λ},\tilde{u})g_f(\mathrm{\λ},\tilde{u})\right]}_{\tilde{u}=\tilde{u}*(\mathrm{\λ},\stackrel{\→}{x})},$

$\tilde{u}*(\mathrm{\λ},\stackrel{\→}{x})=\frac{D\stackrel{\→}{x}.{\stackrel{\→}{e}}_2}{R+\stackrel{\→}{x}.{\stackrel{\→}{e}}_1},{\stackrel{\→}{e}}_1=(-\cos \mathrm{\λ},-\sin \mathrm{\λ}),{\stackrel{\→}{e}}_2=(-\sin \mathrm{\λ},\cos \mathrm{\λ}),$

$w(\mathrm{\λ},\tilde{u})=\frac{c\left(\mathrm{\λ}\right)}{c\left(\mathrm{\λ}\right)+c(\mathrm{\λ}+\mathrm{\π}-2\arctan (\tilde{u}/D))},$

$c\left(\mathrm{\λ}\right)=\left\{\begin{array}{ccc}{\cos }^2\frac{\mathrm{\π}(\mathrm{\λ}-d)}{2d}& \mathrm{if}& 0\≤\mathrm{\λ}\≤d\\ 1& \mathrm{if}& d\≤\mathrm{\λ}\≤{\mathrm{\λ}}_u-d\\ {\cos }^2\frac{\mathrm{\π}(\mathrm{\λ}-{\mathrm{\λ}}_u+d)}{2d}& \mathrm{if}& {\mathrm{\λ}}_u-d\≤\mathrm{\λ}\≤{\mathrm{\λ}}_u\\ 0& \mathrm{if}& \mathrm{others}\end{array}\right.,$

d＝10°

以上所述仅为本发明的较佳实施例而已，并非用来限定本发明的实施范围。任何所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰，因此本发明的保护范围应当视权利要求书所界定范围为准。

Claims (5)

1.一种基于普通C臂机的脊柱横断面图像合成方法，其特征在于：包括以下步骤：

步骤一、对普通C臂机获得的2D锥束投影图像进行预处理，得到1D扇束投影数据；

步骤二、对步骤一得到的扇束进行虚拟旋转再处理，计算得到虚拟扇束投影数据；

步骤三、对2D断层图像合成方法中的相关参数进行调整，采用调整后的方法基于虚拟扇束投影数据合成脊柱横断面图像。

2.如权利要求1所述的基于普通C臂机的脊柱横断面图像合成方法，其特征在于：步骤一中所述的预处理，是指首先对2D锥束投影图像进行中值滤波除噪，然后提取出图像中间的3列或2列像素的灰度值，对提取的像素的灰度值取均值构成1D扇束投影数据；

所述1D扇束投影与普通C臂机获得的2D锥束投影的关系为：扇束光源点与C臂X线锥束光源点重合，扇束中轴线与C臂X线锥束中轴线重合，扇束光源点到其成像直线之间的距离等于C臂X线锥束光源点到其成像平面之间的距离，扇束成像直线位于C臂X线锥束成像平面上，扇束的扇形角等于C臂X线锥束的锥角γ_m。

3.如权利要求1或2所述的基于普通C臂机的脊柱横断面图像合成方法，其特征在于：步骤二中所述的再处理，是指对步骤一所得扇束绕其光源点进行虚拟旋转得到虚拟扇束，通过计算后得到对应的虚拟扇束投影数据；

所述虚拟扇束与扇束的关系为：虚拟扇束光源点到其成像直线之间的距离等于扇束光源点到其成像直线之间的距离D，虚拟扇束光源点到其成像直线的垂线通过C臂旋转中心；对虚拟扇束投影和扇束投影分别建立笛卡尔坐标系，令两个坐标系的原点分别位于光源点在各自成像直线上的垂足，且令虚拟扇束投影坐标轴正向与扇束投影坐标轴正向所夹角度为θ，则两个坐标系下投影像素坐标及灰度值对应关系为：

$u=D.\tan (\arctan \left(\frac{u_v}{D}\right)-\mathrm{\θ}),$ $g_v(\mathrm{\λ},u_v)=g_c(\mathrm{\λ},u).\left(\frac{\sqrt{D^2+u_v^2}}{\sqrt{D^2+u^2}}\right)$

上式中，u为扇束投影在其笛卡尔坐标系中的像素坐标，u_v为虚拟扇束投影在其笛卡尔坐标系中的像素坐标，g_c (λ，u)为扇束投影在对应坐标处的像素灰度值，g_v(λ，u_v)为虚拟扇束投影在对应坐标处的像素灰度值，λ为C臂旋转角度。

4.如权利要求1所述的基于普通C臂机的脊柱横断面图像合成方法，其特征在于：步骤三中所述的对2D断层图像合成方法中的相关参数进行调整，具体包括：将光源点沿圆弧路径运动的半径R调整为光源点到C臂旋转中心的距离，将对扇束投影数据进行滤波处理的区间下限u_l调整为-D.tan(γ_m-θ)，上限u_u调整为D.tan(γ_m+θ)，上述各式中，γ_m为扇束的扇形角，D为扇束光源点到其成像直线之间的距离，θ为虚拟扇束投影坐标轴正向与扇束投影坐标轴正向所夹角度。

5.如权利要求2所述的基于普通C臂机的脊柱横断面图像合成方法，其特征在于：所述提取出图像中间的3列或2列像素的灰度值，是指当图像像素列数n为奇数时，提取出图像中间的3列像素的灰度值，当图像像素列数n为偶数时，提取出图像中间的2列像素的灰度值。

Images