CN102542600B - 一种基于cuda技术的仿真投影drr生成方法 - Google Patents

Abstract

一种基于CUDA技术的仿真投影DRR生成方法。上述方法包括：建立一种新的投影模型；对二维图像的每个像素点进行反投影；通过计算反投影线上三维图像的灰度值之和，最终生成图像配准所需的仿真投影图像。并且在该模型的基础上，本专利还提出了一种基于CUDA(Compute Unified Device Architecture，统一计算设备架构)技术的硬件加速方法。采用本发明所述方法，可以有效地消除投影失真和尺度变换之间的耦合，增强二维与三维图像配准的鲁棒性和精度，并且提高算法的实时性与效率。

Description

一种基于CUDA技术的仿真投影DRR生成方法

技术领域

本发明涉及医学图像配准领域，尤其涉及二维三维医学图像配准中的投影模型以及仿真投影(DRR)的生成方法。

背景技术

计算机辅助的医学手术是一个集医学和计算机技术等诸多学科为一体的新型交叉研究领域。其目的是借助计算机以及跟踪设备模拟手术中涉及的各个过程包括手术规划，手术导航等来指导医生实现高精度的外科手术。

医学图像配准是计算机辅助手术的关键技术。特别是在骨创伤导航手术中，通过二维三维图像配准可以将术前的二维手术规划数据和术中三维图像匹配，进而完成三维图像中手术规划数据的重建，达到指导手术的目的。

二维三维图像配准中，仿真投影(DRR)定义了当前三维图像的位置。通过仿真投影和投影图像比较，可以矫正三维图像的位置最终达到和投影图像的匹配。

图像配准中，传统的投影模型假设投影结构是一个正四面锥。并且X光源到投影平面的距离已知。待优化的参数是投影模型和三维图像的空间位置关系，由三个旋转参数和三个平移参数组成。这种模型对于骨创伤手术中的配准并不适用。首先，当三维图像沿中心投影移动时会同时产生投影失真和尺度变化。这两种变化的耦合会降低配准方法的鲁棒性和精确度。

在生成DRR的过程中，需要进行大量的计算。这影响了配准过程的实时性。通过以CUDA为平台的GPU硬件运算可以增强DRR生成的实时性。

发明内容

为了克服由于投影模型带来的配准鲁棒性差和精度低的缺点。本发明提出了一种新的投影模型以及基于这个模型的仿真投影生成方法和基于CUDA平台的实现。

根据本发明的一个方面，提出了应用于二维三维图像配准的投影模型建立方法包括以下步骤：

A01：以三维图像的中心为坐标原点建立笛卡尔坐标系即三维图像坐标系XYZ，使得X，Y，Z轴分别与三维图像中对应的外平面正交，记坐标原点为O；

A02：以三维图像坐标系XYZ的原点为坐标原点建立笛卡尔坐标系即投影模型坐标系UVW，使得坐标轴U，V，W与坐标轴X，Y，Z分别保持一致，记坐标原点为ISO；

A03：为投影模型建立一个X光源，X光源位于W轴正半轴上，距投影模型坐标系UVW坐标原点ISO距离为D1；

A04：在接收器所在的平面为投影模型建立一个投影面；投影面正交于W轴；W轴与投影面的交点位于投影面几何中心，且位于W轴负半轴；交点与投影模型坐标系UVW的原点ISO距离为D2；

A05：为投影模型建立一个虚拟投影面；虚拟投影面平行于投影面且中心通过投影模型坐标系UVW的坐标原点ISO；

A06：将投影模型相对三维图像的运动描述为投影模型坐标系UVW在三维图像坐标系XYZ内的旋转和平移，包括：平面外绕U，V轴的旋转Ru、Rv，平面内绕W轴的旋转参数Rw，平面内平移Tu、Tv即UVW在XYZ坐标系内沿UV平面的平移；这个运动记为^XYZT_UVW； $T_{\mathrm{UVW}}^{\mathrm{XYZ}}=\left[\begin{array}{cccc}\cos \mathrm{\β}& 0& \sin \mathrm{\β}& 0\\ 0& 1& 0& 0\\ -\sin \mathrm{\β}& 0& \cos \mathrm{\β}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \mathrm{\α}& -\sin \mathrm{\α}& 0\\ 0& \sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}\cos \mathrm{\φ}& -\sin \mathrm{\φ}& 0& 0\\ \sin \mathrm{\φ}& \cos \mathrm{\φ}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& \mathrm{Tu}\\ 0& 1& 0& \mathrm{Tv}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$ α为投影模型坐标系UVW围绕U轴旋转角度Ru，β为投影模型坐标系UVW围绕V轴旋转角度Rv，φ为投影模型坐标系UVW围绕W轴旋转角度Rw。>

A07：投影模型自身的变化由X光源到三维图像距离D1和三维图像到接收器距离D2描述，记为I(D1，D2)；

根据本发明的另一个方面，DRR的生成方法包括以下步骤：

B01：对于投影平面上的任意点的坐标记为^UVWP_D(u_D,v_D,-D₂)，其中u_D为虚拟平面上任意点的U轴坐标，v_D为投影平面上任意点的v轴坐标，D₂为UVW坐标原点ISO到接收器距离；

对于投影面上的任意点^UVWP_D(u_D,v_D,-D₂)求其在虚拟投影面上的对应点的坐标为^UVWP(u,v,0)其中 $\left\{\begin{array}{c}u=u_D\frac{D_1}{D_1+D_2}\\ v=v_D\frac{D_1}{D_1+D_2}\end{array}\right.;$ D₁是X光源到UVW坐标原点ISO的距离；D₂是UVW坐标原点ISO到接收器的距离；u_D是投影面上的点在U坐标轴的坐标值；v_D是投影面上的点在V坐标轴的坐标值；u是虚拟投影面上的对应点在U坐标轴的坐标值；v是虚拟投影面上的对应点在V坐标轴的坐标值；

B02：对由坐标为^UVWP(u,v,0)的点和坐标为^UVWP(0,0,D₁)的X光源所在点所确定的直线上的所有点上灰度值求和得到虚拟投影面上点^UVWP(u,v,0)的灰度值；这条直线上任意点可以表示为其中^UVWP_l表示该直线上的点在UVW坐标系中的坐标值；w表示直线上点在W轴上的坐标；u表示该直线与虚拟投影面交点在U轴上的坐标；v表示该直线与虚拟投影面交点在V轴上的坐标；D₁是X光源到UVW坐标原点ISO的距离；

B03：将坐标^UVWP_l变换到XYZ坐标系下得到该点坐标为^XYZP_l＝^XYZT_UVW*^UVWP_l；其中^XYZT_UVW是坐标系XYZ到坐标系UVW的转换矩阵；^UVWP_l是直线上点在坐标系UVW下的坐标；^XYZP_l是直线上的点在坐标系XYZ下的坐标；将坐标值^XYZP_l四舍五入取整得^XYZP_{l_Z}；设三维图像的密度函数为G(x,y,z)，则二维投影上点^UVWP_D(u_D,v_D,-D₂)的灰度值可以由下面的等式得出；

$I(u_D,v_D,-D_2)=\underset{w=-\frac{s}{2}}{\overset{\frac{s}{2}}{\mathrm{\Σ}}}G\left(P_{l\_Z}^{\mathrm{XYZ}}\right)$

其中s是三维图像对角线的长度；I(u_D,v_D,-D₂)是沿着当前射线的总吸收量；u_D表示该直线与投影面交点在U轴上的坐标；v_D表示该直线与虚拟投影面交点在V轴上的坐标；D₂是ISO中心到接收器的距离，^XYZP_{l_Z}是直线上点在XYZ坐标系下取整后的坐标；

B04：根据Beer’s定理投影平面上对应点的灰度值其中H_max是投影图像的最大灰度值；I(u_D,v_D,-D₂)是沿着当前射线的总吸收量，这里的e是数学常数，就是自然对数的底数，近似等于2.718281828；

B05：对接收器上每一个像素重复步骤B01-B04得到仿真投影DRR；

根据本发明的另一方面，DRR生成方法由以下步骤：

初始化CUDA；

分配内存；

三维图像传输到显卡设备；

初始化MxN个线程，其中M和N代表输出图像维度。每一个线程对应着图像上的一个像素；

按照B-E中描述的方法求当前像素灰度值；

投影图像输出到主程序。

本发明可以获得如下有益效果：

1.本方法的投影模型通过分离X光源-ISO距离D1和ISO-接收器距离D2有效地去除投影失真和尺度优化之间的耦合，从而提高配准的鲁棒性和精确度。

2.本方法对仿真投影的生成进行硬件加速，可以实时生成仿真投影图像

3、本方法提出适用于骨科创伤手术图像配准的投影模型。

附图说明

图1示出投影模型以及参数。

图2示出基于CUDA技术的DRR生成流程

图3示出投影面上每一像素灰度值的确定方法

具体实施方式

下面结合附图和具体实例对本发明提出的一种基于CUDA平台的仿真投影图像生成方法进行详细描述。

如图一所示，投影模型由三维图像和投影结构两部分组成，他们分别在XYZ直角坐标系和UVW直角坐标系下被定义。其中三维图像的中心位于直角坐标系的原点O。投影结构以坐标系UVW为参考坐标系。设UVW的原点为ISO中心。X光源位于W轴正半轴上。X光源到ISO中心的距离是D1。投影平面垂直于W轴。投影平面和W轴的交点位于W轴的负半轴。交点位于投影平面的中心。交点和ISO中心的距离是D2。定义虚拟平面是过ISO中心并且和投影平面平行的平面。W轴与虚拟平面的交点位于虚拟平面的中心。

UVW在坐标系XYZ下的空间变化可以由以下变化组成：

围绕V轴的旋转Rv→围绕U轴的旋转Ru→围绕W轴的旋转Rw→沿着UV平面的平移Tu，Tv。

其中围绕U，V轴的旋转为平面外旋转，影响着投影中物体的形状；围绕W轴的旋转为平面内旋转，影响着投影中物体的方向；Tu和Tv是平面内的平移，影响着投影中物体的位置。^XYZT_UVW为坐标系XYZ到坐标系UVW的转换矩阵。

$T_{\mathrm{UVW}}^{\mathrm{XYZ}}=\left[\begin{array}{cccc}\cos \mathrm{\β}& 0& \sin \mathrm{\β}& 0\\ 0& 1& 0& 0\\ -\sin \mathrm{\β}& 0& \cos \mathrm{\β}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \mathrm{\α}& -\sin \mathrm{\α}& 0\\ 0& \sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}\cos \mathrm{\φ}& -\sin \mathrm{\φ}& 0& 0\\ \sin \mathrm{\φ}& \cos \mathrm{\φ}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& \mathrm{Tu}\\ 0& 1& 0& \mathrm{Tv}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$

α为投影模型坐标系UVW围绕U轴旋转角度Ru，β为投影模型坐标系UVW围绕V轴旋转角度Rv，φ为投影模型坐标系UVW围绕W轴旋转角度Rw。

投影结构的自身变化可以由参数D1，D2来决定。D1的变化影响着投影中物体的投影失真；D2的变化影响着投影中物体的尺度。投影平面上的点^UVWP_D(u_D,v_D,-D₂)和虚拟投影面上的点^UVWP(u,v,0)存在着如下关系:

$\left\{\begin{array}{c}u=u_D\frac{D_1}{D_1+D_2}\\ v=v_D\frac{D_1}{D_1+D_2}\end{array}\right.$

由点^UVWP(u,v,0)以及X光源可以确定一条X光射线。射线上的点可以描述为其中u，v是射线和虚拟投影面交点的坐标。w是射线上的点在W轴方向上的坐标值，D1为X光源到UVW坐标原点ISO的距离，D2为UVW坐标原点到接收器的距离。

如图二所示，首先初始化CUDA平台。接着分配处理中所需的内存包括：一个3D Texture用来存放三维图像；一个2D Texture用来存放生成的二维投影图像。在显卡设备中分配MxN个处理单元。其中M，N是投影图像的维度。这样每一个处理单元对应着图像中的一个像素。接下来对每一个像素按照图三所示的方法求灰度值。如图三示出，投影平面上点^UVWP_D(u_D,v_D,-D₂)的强度可以通过对射线上所有点的强度求和得到。考虑到效率这里只对以w＝0为中心区间[-s/2,s/2]之间的强度进行求和。其中s是三维图像对角线的长度。为了求对角线上点的灰度值，先将UVW坐标系下点坐标通过公式^XYZP_l＝^XYZT_UVW*^UVWP_l变换到XYZ坐标系下。再将点^XYZP_l各个分量的坐标以四舍五入的方式取整得到^XYZP_{l_Z}。设三维图像在坐标X，Y，Z点的灰度值是G(X,Y,Z)。则取整后的点^XYZP_{l_Z}的灰度值为G(^XYZP_{l_Z})。接下来对范围 $w\∈\left[\begin{array}{cc}-\frac{s}{2}& \frac{s}{2}\end{array}\right]$ 内射线上所有强度值求和得到投影面上点^UVWP_D(u_D,v_D,-D₂)的强度值其中s是三维图像对角线的长度，^XYZP_{l_Z}是直线上点在XYZ坐标系下取整后的坐标。最后，用Beer’s理论求出该像素的灰度值其中H_max是投影图像的最大灰度值；I(u_D,v_D,-D₂)是沿着当前射线的总吸收量，这里的e是数学常数，就是自然对数的底数，近似等于2.718281828。求出的H(u_D,v_D,-D₂)就是投影平面上点^UVWP_D(u_D,v_D,-D₂)的灰度值。在求出所有平面上点灰度值之后，将投影图像传回主程序完成DRR的生成。

Claims (2)

1.一种基于CUDA技术的仿真投影DRR生成方法,其特征在于：该方法包括如下步骤：

A、建立投影模型，包括如下步骤：

A01：以三维图像的中心为坐标原点建立笛卡尔坐标系即三维图像坐标系XYZ，使得X，Y，Z轴分别与三维图像中对应的外平面正交，记坐标原点为O；

A02：以三维图像坐标系XYZ的原点为坐标原点建立笛卡尔坐标系即投影模型坐标系UVW，使得坐标轴U，V，W与坐标轴X，Y，Z的方向分别保持一致，记坐标原点为ISO；

A03：为投影模型建立一个X光源，X光源位于W轴正半轴上，距投影模型坐标系UVW坐标原点ISO距离为D1；

A04：为投影模型建立一个投影面；投影面与接收器共面、大小相等且方向相同；投影面正交于W轴；W轴与投影面的交点位于投影面几何中心，且位于W轴负半轴；交点与投影模型坐标系UVW的原点ISO距离即投影模型坐标系UVW坐标原点ISO到投影面距离为D2；投影模型坐标系UVW的U轴和V轴分别与投影面的对应边平行；

A05：为投影模型建立一个虚拟投影面；虚拟投影面平行于投影面且中心通过投影模型坐标系UVW的坐标原点ISO；

A06：将投影模型相对三维图像的运动描述为投影模型坐标系UVW在三维图像坐标系XYZ内的旋转和平移，包括：平面外绕U，V轴的旋转Ru、Rv，平面内绕W轴的旋转参数Rw，平面内平移Tu、Tv即UVW在XYZ坐标系内沿UV平面的平移；这个运动记为^XYZT_UVW；

$T_{\mathrm{UVW}}^{\mathrm{XYZ}}=\left[\begin{array}{cccc}\cos \mathrm{\β}& 0& \sin \mathrm{\β}& 0\\ 0& 1& 0& 0\\ -\sin \mathrm{\β}& 0& \cos \mathrm{\β}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \mathrm{\α}& -\sin \mathrm{\α}& 0\\ 0& \sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}\cos \mathrm{\φ}& -\sin \mathrm{\φ}& 0& 0\\ \sin \mathrm{\φ}& \cos \mathrm{\φ}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& \mathrm{Tu}\\ 0& 1& 0& \mathrm{Tv}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$ α为投影模型坐标系UVW围绕U轴旋转角度Ru，β为投影模型坐标系UVW围绕V轴旋转角度Rv，φ为投影模型坐标系UVW围绕W轴旋转角度Rw；

A07：投影模型自身的变化由X光源到投影模型坐标系UVW坐标原点ISO距离D1和投影模型坐标系UWV坐标原点ISO到投影面距离D2描述，记为I（D1，D2）；

B、DRR的生成方法包含如下步骤：

B01：对于投影面上的任意点的坐标记为^UVWP_D(u_D,v_D,-D₂)，其中u_D为投影面上的点的U轴坐标，v_D为投影面上的点的V轴坐标，D₂为UVW坐标原点ISO到投影面距离；

对于投影面上的任意点^UVWP_D(u_D,v_D,-D₂)求其在虚拟投影面上的对应点的坐标为^UVWP(u,v,0)其中 $\left\{\begin{array}{c}u=u_D\frac{D_1}{D_1+D_2}\\ v=v_D\frac{D_1}{D_1+D_2}\end{array}\right.;$ D₁是X光源到UVW坐标原点ISO的距离；D₂是UVW坐标原点ISO到投影面的距离；u_D是投影面上的点在U坐标轴的坐标；v_D是投影面上的点在V坐标轴的坐标；u是虚拟投影面上的点在U坐标轴的坐标；v是虚拟投影面上的点在V坐标轴的坐标；

B02：对由坐标为^UVWP(u,v,0)的点和坐标为^UVWP(0,0,D₁)的X光源所在点所确定的直线上的所有点上灰度值求和得到虚拟投影面上点^UVWP(u,v,0)的灰度值；这条直线上任意点可以表示为其中^UVWP_l表示该直线上的点在UVW坐标系中的坐标值；w表示直线上点在W轴上的坐标；u表示虚拟投影面上的点在U轴上的坐标；v表示虚拟投影面上的点在V轴上的坐标；D₁是X光源到UVW坐标原点ISO的距离；

B03：将坐标^UVWP_l变换到XYZ坐标系下得到该点坐标为^XYZP_l=^XYZT_UVW*^UVWP_l；其中^XYZT_UVW是坐标系XYZ到坐标系UVW的转换矩阵；^UVWP_l是直线上点在坐标系UVW下的坐标；^XYZP_l是直线上的点在坐标系XYZ下的坐标；将坐标值^XYZP_l四舍五入取整得^XYZP_{l_Z}；设三维图像的密度函数为G(x,y,z)，则二维投影上点^UVWP_D(u_D,v_D,-D₂)的灰度值可以由下面的等式得出；

$I(u_D,v_D,-D_2)=\underset{w=-\frac{s}{2}}{\overset{\frac{s}{2}}{\mathrm{\Σ}}}G\left(P_{l\_Z}^{\mathrm{XYZ}}\right)$

其中其中s是三维图像对角线的长度；I(u_D,v_D,-D₂)是沿着当前射线的总吸收量；u_D表示投影面上的点在U轴上的坐标；v_D表示投影面上的点在V轴上的坐标；D₂是UVW坐标原点ISO到投影面的距离，^XYZP_{l_Z}是直线上点在XYZ坐标系下取整后的坐标；

B04：根据Beer’s定理投影面上对应点的灰度值其中H_max是投影图像的最大灰度值；I(u_D,v_D,-D₂)是沿着当前射线的总吸收量，e是自然对数的底数；

B05：对投影面上每一个像素重复步骤B01-B04得到仿真投影DRR。

2.如权利要求1所述的DRR生成方法的流程，其特征在于：其步骤有：

1）初始化CUDA；

2）分配内存；

3）传输三维图像到显卡设备；

4）在显卡设备中初始化MxN个线程，其中M和N代表输出图像维度；每一个线程对应着图像上的一个像素；

5）对投影面上每一像素求灰度值；

6）输出投影图像到主程序。