CN102945328A - 基于gpu并行运算的x射线造影图像仿真方法 - Google Patents

Abstract

本发明提出了一种基于GPU并行运算的X射线造影图像仿真方法，能够提高仿真结果的真实度及仿真算法的时间效率。主要包括1、构建旋转造影成像模型；2、利用区域增长的方法对CT影像的冠脉进行分割，通过冠脉CT值的权重优化，模拟冠脉增强的三维影像；3、根据旋转造影成像系统的参数设定，建立透视投影模型，计算三维影像空间到二维成像空间的映射关系；4、构建能量衰减模型，计算透射X射线的能量分布；5、以三维模拟影像的体素坐标对GPU的纹理数据初始化，利用透视投影模型和能量衰减模型并行计算仿真造影图像。

Description

基于GPU并行运算的X射线造影图像仿真方法

技术领域

本发明涉及一种基于GPU并行运算的X射线造影图像仿真方法，适用于重建算法优势的评价及临床医务人员培训等领域。

背景技术

冠状动脉造影成像以其高分辨率和高对比度等特征，被称为心血管疾病诊断和治疗成像方法中的“金标准”。获取造影图像的途径包括单平面、双平面及旋转造影成像方式。旋转造影成像可实现多角度照射，降低了对狭窄病灶误诊的概率，方便医生对血管狭窄程度的了解及合适治疗方案的制定。但临床冠脉造影图像仍普遍存在以下问题：(1)因血管结构错综复杂引起的交叉和遮挡严重；(2)二维造影图像丢失大量的空间信息；(3)造影图像中存在大量的背景噪声。上述问题导致临床医生很难通过有限的成像角度对病灶进行诊断和治疗。

基于X射线冠状动脉造影图像的三维重建可有效弥补冠状动脉造影术中存在的缺陷。该技术可从有限的二维图像中恢复血管的三维信息，在先验知识的引导下，实现对血管狭窄部分的测量和评估、对冠状动脉运动的模拟、及最佳造影成像角度的确定等。但因真实造影数据的多样性导致三维重建算法缺少统一的评价标准。

基于X射线造影图像的仿真可获取具有先验知识的造影图像，为三维重建算法的研究和改进提供精确的评价标准，最终实现基于真实造影图像高精度的三维重建。在临床应用中，有利于医生完成对病灶的快速定位和处理，有效地减少手术时间，提高手术效率，减少X射线的辐射对病人和医生造成的危害。

实现X射线造影图像仿真的主流算法可分为两类——体素驱动法和光线驱动法。体素驱动法是通过连接三维影像的所有体素和光源，利用透视投影模型，获取体素在投影平面上对应的像素坐标及所有像素的灰度分布信息，从而实现造影图像的仿真，但该方法因精度较低而未被广泛使用。光线驱动法是通过连接投影平面上的像素和光源，利用能量积分衰减模型，求解像素的灰度值，获取仿真的造影图像。该方法已被广泛用于图像去噪及重建的后续处理中，但其计算量偏大，故在实际算法设计中，可根据精度要求，设定固定的能量积分步长，获取造影图像，因而能在一定程度上降低计算的复杂度。随着基于GPU的单指令多数据计算模型的迅速发展，在医学影像的处理中，通过设计并行算法程序，以体素或像素初始化GPU纹理，实现并行渲染，可极大程度的提高算法的运行效率。

综上所述，现有X射线造影图像仿真存在的问题：

1.不能对造影图像进行全方位成像，实现病患的全面性诊断；

2.仿真结果和真实造影图像相比，真实度低；

3.仿真造影图像时间效率低。

发明内容

为克服现有X射线造影图像仿真技术的不足，本发明提出了一种基于GPU并行运算的X射线造影图像仿真方法，是一种基于冠脉增强的三维影像，以透视投影模型和能量衰减模型作为成像基础的造影图像的快速仿真方法。

该种基于GPU并行运算的X射线造影图像仿真方法，包括以下步骤：

第一步：旋转造影成像模型的构建：根据临床C型臂旋转造影成像系统的参数信息，包括C型臂旋转角度、成像系统各部位之间的距离，构建旋转造影成像模型；

所述的旋转造影成像支持静态和动态两种成像模式，其中静态模式通过固定成像角度，实现病患的实时成像；动态模式通过设置成像旋转范围，实现病患的全方位成像；

第二步：三维影像中冠状动脉的分割：通过选取血管种子点，利用区域增长算法对CT影像中的冠状动脉进行半自动分割；

第三步：基于权重优化的冠状动脉增强影像的模拟：根据冠脉结构在三维影像中的空间分布信息，对相应的CT值进行权重优化，增强三维影像中冠状动脉脉络结构；

第四步：透视投影模型的构建：根据摄像机成像原理，利用旋转成像系统的参数信息，构建内参和外参矩阵，实现三维影像空间到二维成像空间的映射；

第五步：基于欧式距离的权重因子设定：利用第四步获取的三维影像的投影坐标信息，计算每个投影坐标点和其最近四个像素坐标之间的欧氏距离，构建归一化比例因子，作为不同体素对成像平面像素灰度信息的权值；

第六步：能量透射衰减模型的构建：根据X射线在三维影像中透射路径信息，构建指数型能量衰减积分模型，计算透射X射线的能量分布；

第七步：GPU纹理参数的初始化：将三维影像的每层体素坐标设置为GPU加速的相关纹理参数，作为后续算法并行处理的原始数据；

第八步：基于并行计算的冠脉增强的X射线造影成像仿真：利用第七步所设定的纹理参数，采用第四步、第五步和第六步所述的X射线造影成像模型，并行计算三维影像在投影平面的灰度信息分布，实现最终冠脉增强的造影仿真影像。

本发明的有益效果：

与现有方法相比，本方法的优点在于利用权重优化的冠脉造影图像的模拟，可提高冠脉的对比度，且基于GPU的并行计算大大提高算法的运行效率。

附图说明

图1是本发明X射线造影图像仿真方法流程图；

图2是X射线旋转冠脉造影模拟系统图。

具体实施方式

下面结合具体实施例和附图详细说明本发明，但本发明并不仅限于此。

图1为本发明的流程图，所述的基于GPU并行计算的造影图像的仿真包括以下几个步骤：

步骤S101，根据临床C型臂旋转造影成像系统的参数信息，包括C型臂旋转角度、成像系统各部位之间的距离，构建旋转造影成像模型，如图2所示。具体内容如下：

临床上基于C型臂的旋转造影成像系统被用于获取X射线造影图像。在该系统中，X射线源和成像平面分别固定于C型臂的两端，C型臂的旋转中心和病患的中心一致。旋转造影系统支持静态和动态两种成像模式，其中静态模式通过固定成像角度，实现病患的实时成像，常用于心脏介入手术；动态模式则通过设置成像旋转范围，实现病患的全方位成像，可用于病患的诊断。

在旋转造影成像系统中，C型臂绕旋转中心P带动光源和成像平面旋转。射线源到旋转中心和成像平面的距离分别为SOD和SID，SOD描述了成像系统的旋转半径，参数标注如图2所示，θ和

为C型臂旋转的主角度和次角度，其中θ描述成像时的旋转角度，代表了左前斜角(RAO)和右前斜角(LAO)，旋转范围为RAO120°~LAO180°，

描述射线源近头(CRA)或近脚(CAU)的位置，变化范围为CRA25°~CAU25°。在成像时，影像的中心和成像系统的旋转中心始终保持一致，以保证投影图像始终位于成像平面中。影像的大小为成像平面的大小为

步骤S102，通过选取血管种子点，利用区域增长算法对CT影像中的冠状动脉进行半自动分割。

步骤S103，根据冠脉结构在三维影像中的空间分布信息，对相应的CT值进行权重优化，增强三维影像中冠状动脉脉络结构。如下式所述：

步骤S104，根据摄像机成像原理，利用旋转成像系统的参数信息，构建内参和外参矩阵，形成透射投影模型，实现三维影像空间到二维成像空间的映射。具体过程如下：

成像系统中包含三个坐标系：世界坐标系o-x_Wy_Wz_W，投影图像坐标系o’-uv和摄像机坐标系c-x_Cy_Cz_C，其中世界坐标系以影像中心为原点，以实际物理距离为单位；投影图像坐标系以投影平面起始位置为原点，以像素为单位；摄像机坐标系以光源中心为原点，以体素为单位。具体坐标系分布如图2示。

三个坐标系中的齐次坐标分别表示如下：

世界坐标系：p_W＝[x_W,y_W,z_W,1]^T

摄像机坐标系：p_C＝[x_C,y_C,z_C,1]^T

投影图像坐标系：p_I＝[u,v,1]^T

通过左乘变换矩阵可实现各坐标系间的相互变换。设M_C，M_I表示坐标系间的变换矩阵，且分别实现o-x_Wy_Wz_W坐标系到c-x_Cy_Cz_C坐标系的变换，和c-x_Cy_Cz_C坐标系到o’-uv坐标系的变换：

p_c＝M_cp_W

p_I＝M_Ip_c

则有p_I＝M_IM_Cp_W，在摄像机成像模型中对应有： $p_I=K\left[\begin{array}{cc}R& t\\ 0_3^T& 1\end{array}\right]p_W/z_C.$ 其中K是摄像机的内参矩阵， $\left[\begin{array}{cc}R& t\\ 0_3^T& 1\end{array}\right]$ 是摄像机的外参矩阵。p_W表示在世界坐标系中，体素的原始坐标相对X射线源的坐标点对称处理后的坐标，z_C为体素在摄像机坐标系中的Z坐标。

在摄像机成像模型中，矩阵K的表达形式如下：

$K=\left[\begin{array}{cccc}f/\mathrm{dx}& 0& u_0& 0\\ 0& f/\mathrm{dy}& v_0& 0\\ 0& 0& 1& 0\end{array}\right]$

在旋转造影仿真系统中，f对应光源到成像平面的距离，即SID；dx和dy为三维影像在平行XOY平面方位的分辨率。

矩阵R和t分别描述世界坐标系到摄像机坐标系的旋转和平移。表示如下：

其中θ为C型臂与XOZ平面的夹角，

为C型臂与YOZ平面的夹角，SOD为射线源到旋转中心P的距离。此时公式描述中设定X射线源在世界坐标系的第四象限，在其它象限只需重新设定对应的平移矩阵t。

投影模型可表示如下：

步骤S105，利用步骤S104获取的三维影像的投影坐标信息，计算每个投影坐标点和其最近四个像素坐标之间的欧氏距离，构建归一化比例因子，作为不同体素对成像平面像素灰度信息的权值。权值计算公式如下：

${\mathrm{\λ}}_i=\underset{k}{\min }\frac{d_k}{(d1+d2+d3+d4)}(k=\mathrm{1,2,3,4})$

其中λ_i的倒数表示第i个体素对当前像素灰度值的贡献，d₁,d₂,d₃,d₄分别表示每个投影坐标点和其最近四个像素坐标之间的欧氏距离。

步骤S106，根据X射线在三维影像中透射路径信息，构建指数型能量衰减积分模型，计算透射X射线的能量分布。具体方法如下：

X射线穿过三维影像发生能量衰减，能量衰减计算函数表示如下：

I＝I₀e^-∫udi

其中I₀和I分别表示光线的初始能量和衰减后的能量，μ为透射路径上的衰减系数，计算如下：

$\mathrm{\μ}=[\frac{\mathrm{CT}}{1000}+1{]\mathrm{\μ}}_w$

其中CT表示当前体素的CT值，单位为HU(Hounsfield Unit)，人体不同组织的CT值处于不同的范围；μ_w表示水的衰减系数，其值随入射X射线能量的增加而减小。

根据步骤S105计算的贡献权值分布，累加满足条件的衰减后的能量值，即可获得像素的灰度值：

下标i表示对像素P灰度值有贡献的体素的编号。r^Z是三维影像沿Z轴方向的分辨率，

对应光线在体素内的途经长度。

步骤S107，将三维影像的每层体素坐标设置为GPU加速的相关纹理参数，作为后续算法并行处理的原始数据。

在该步骤中，采用OpenGL Extension中提供的矩形纹理(Texture Rectangles)模式，允许数据以任意大小存储为纹理。

步骤S108，利用步骤S107所设定的纹理参数，采用步骤S104、S105和S106所述的X射线造影成像模型，并行计算三维影像在投影平面的灰度信息分布，实现最终冠脉增强的造影仿真影像。

将利用步骤S105和S106所得纹理的输出结果作为中间数据，存储在纹理渲染(Render toTexture)的缓存中，作为步骤S107的输入数据，计算投影平面的灰度信息分布，则获得最终的冠脉增强的造影仿真影像。

虽然参考优选实施例对本发明进行描述，但以上所述实例并不构成本发明保护范围的限定，任何在本发明的精神及原则内的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围内。

Claims (4)

1.基于GPU并行运算的X射线造影图像仿真方法，其特征在于，包括以下步骤：

第一步：旋转造影成像模型的构建：根据临床C型臂旋转造影成像系统的参数信息，包括C型臂旋转角度、成像系统各部位之间的距离，构建旋转造影成像模型；

第二步：三维影像中冠状动脉的分割：通过选取血管种子点，利用区域增长算法对CT影像中的冠状动脉进行半自动分割；

第三步：基于权重优化的冠状动脉增强影像的模拟：根据冠脉结构在三维影像中的空间分布信息，对相应的CT值进行权重优化，增强三维影像中冠状动脉脉络结构；

第四步：透视投影模型的构建：根据摄像机成像原理，利用旋转成像系统的参数信息，构建内参和外参矩阵，实现三维影像空间到二维成像空间的映射；

第五步：基于欧式距离的权重因子设定：利用第四步获取的三维影像的投影坐标信息，计算每个投影坐标点和其最近四个像素坐标之间的欧氏距离，构建归一化比例因子，作为不同体素对成像平面像素灰度信息的权值；

第六步：能量透射衰减模型的构建：根据X射线在三维影像中透射路径信息，构建指数型能量衰减积分模型，计算透射X射线的能量分布；

第七步：GPU纹理参数的初始化：将三维影像的每层体素坐标设置为GPU加速的相关纹理参数，作为后续算法并行处理的原始数据；

第八步：基于GPU加速的X射线造影成像：利用第七步所设定的纹理参数，采用第四步、第五步和第六步所述的X射线造影成像模型，并行计算三维影像在投影平面的灰度信息分布，实现最终冠脉增强的造影仿真影像。

2.如权利要求1所述的基于GPU并行运算的X射线造影图像仿真方法，其特征在于，该旋转造影成像系统支持静态和动态两种成像模式，其中静态模式通过固定成像角度，可实现病患的实时成像；动态模式通过设置成像旋转范围，可实现病患的全方位成像。

3.如权利要求1或2所述的基于GPU并行运算的X射线造影图像仿真方法，其特征在于，所述的透视投影模型采用以下方法构建：首先模拟成像系统中的三个坐标系：世界坐标系o-x_Wy_Wz_W，投影图像坐标系o′-uv和摄像机坐标系c-x_Cy_Cz_C，其中世界坐标系以影像中心为原点，以实际物理距离为单位；投影图像坐标系以投影平面起始位置为原点，以像素为单位；摄像机坐标系以光源中心为原点，以体素为单位；

三个坐标系中的齐次坐标表示如下：世界坐标系：p_W＝[x_W，y_W，z_W，l]^T，摄像机坐标系：p_C＝[x_C，y_C，z_C，l]^T，投影图像坐标系：p_I＝[u，v，l]^T；通过左乘变换矩阵实现各坐标系间的相互变换；设M_C，M_I表示坐标系间的变换矩阵，且分别实现o-x_Wy_Wz_W坐标系到c-x_Cy_Cz_C坐标系的变换，和c-x_Cy_Cz_C坐标系到o’-uv坐标系的变换：p_C＝M_Cp_W    p_I＝M_Ip_C；则有P_I＝M_IM_Cp_W，在摄像机成像模型中相应有： $p_I=K\left[\begin{array}{cc}R& t\\ 0_3^T& 1\end{array}\right]p_W/z_C,$ 其中K是摄像机的内参矩阵， $\left[\begin{array}{cc}R& t\\ 0_3^T& 1\end{array}\right]$ 是摄像机的外参矩阵；p_W表示在世界坐标系中，体素的原始坐标相对X射线源的坐标点对称处理后的坐标，z_C为体素在摄像机坐标系中的Z坐标；

在摄像机成像模型中，矩阵K的表达形式如下：

$K=\left[\begin{array}{cccc}f/\mathrm{dx}& 0& u_0& 0\\ 0& f/\mathrm{dy}& v_0& 0\\ 0& 0& 1& 0\end{array}\right]$

在旋转造影仿真系统中，f对应光源到成像平面的距离，即SID；dx和dy为三维影像在平行XOY平面方位的分辨率；

矩阵R和t分别描述世界坐标系到摄像机坐标系的旋转和平移，表示如下：

其中θ为C型臂与XOZ平面的夹角，

为C型臂与YOZ平面的夹角，SOD为射线源到旋转中心P的距离，此时公式描述中设定X射线源在世界坐标系的第四象限，在其它象限只需重新设定对应的平移矩阵t；

投影模型表示如下：

4.如权利要求1或2所述的基于GPU并行运算的X射线造影图像仿真方法，其特征在于，对三维影像中冠脉进行增强，使得后续造影图像模拟中冠脉的对比度提高，可提高临床对血管病灶的确诊率。