CN102663803B - 一种基于RayCasting改进算法的仿真投影DRR生成方法 - Google Patents

Abstract

一种基于RayCasting改进算法的仿真投影DRR生成方法，属医学图像处理领域。该方法是利用三维图像中局部区域的均匀性来判别积分的步长。并通过多尺度图像处理对RayCasting算法实现进一步的加速。本发明在对数据进行预处理时生成用来保存三维数据内容均匀性程度的数据结构，在实际计算光路的过程中根据图像内容的均匀性实现采样点的跳跃，大大提高了DRR生成的速度。

Description

一种基于RayCasting改进算法的仿真投影DRR生成方法

技术领域

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

背景技术

计算机辅助的医学手术导航系统是一个集医学和计算机技术等诸多学科为一体的新型交叉研究领域。其目的是借助计算机以及跟踪设备模拟手术中涉及的各个过程，包括手术规划，手术导航等来指导医生实现高精度的外科手术。随着手术导航系统在临床中的广泛应用，二维三维图像配准已经成为系统开发中的一个关键技术。在二维三维医学图像配准过程中，DRR作为一种常用的中间模型，其生成速度和质量往往对配准的效率和精确度有至关重要的影响(DRR：digitally reconstructed radiography数字重建图像)。

Raycasting算法作为一种典型的像序算法，凭借其较高的图像质量在生成DRR的实际应用中被广泛采用。然而Raycasting算法在DRR计算过程中对像素的遍历以及对光路上像素的积分需要消耗大量的时间，这会大大地降低术中二维三维图像配准的实时性。因此增加算法效率，成为Raycasting算法改进的一个焦点。本发明提出一种基于图像内容均匀性的Raycasting改进算法，在对数据进行预处理时生成用来保存三维数据内容均匀性程度的数据结构，在实际计算光路的过程中根据图像内容的均匀性实现采样点的跳跃，大大改善了DRR生成的速度。

发明内容

本发明提出了一种基于RayCasting改进算法的DRR生成方法。方法包括以下步骤：

A将三维图像V由DICOM服务器调入工作PC，三维图像V的大小是Width×Height×Depth，其中Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，

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

B01：以三维图像的中心为坐标原点建立笛卡尔坐标系即三维图像坐标系XYZ，使得X，Y，Z轴分别与三维图像中对应的外平面正交，记坐标原点为O；三维图像V中像素的坐标由(x，y，z)表示；

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

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

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

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

B06：将投影模型相对三维图像的运动描述为投影模型坐标系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；

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

CDRR的生成方法包含如下步骤：

C01：计算三维图像V的缩小图像V₁...V_n...V_N，其中V_n代表第n级的缩小三维图像，其大小为1≤n≤N，n是当前缩小级数，N是缩小级数的最大值，第n级缩小三维图像V_n中像素的坐标由(x_n，y_n，z_n)表示，第n级的缩小三维图像V_n中像素和三维图像V中的像素存在对应关系V_n(x_n，y_n，z_n)＝V(2ⁿ·x_n，2ⁿ·y_n，2ⁿ·z_n)；

C02：生成每一个缩小三维图像V_n的均匀性矩阵V_{Homo_n}，其大小为1≤n≤N，N是缩小级数的最大值，第n级缩小三维图像V_n中坐标为(x_n，y_n，z_n)的像素的均匀性值为 $V_{\mathrm{Homo}\_n}(x_n,y_n,z_n)=\underset{i=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{j=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{k=-\δ}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\frac{{(V_n(x_n+i,y_n+j,z_n+k)-\stackrel{\‾}{V_{x_n{y}_n{z}_n\mathrm{\δ}}})}^2}{{(2\·\mathrm{\δ}+1)}^3},$ 其中δ是半窗长，是第n级缩小三维图像V_n中以像素(x_n，y_n，z_n)为中心，δ为半径覆盖区域像素的平均值 $\stackrel{\‾}{V_{x_n{y}_n{z}_n\mathrm{\δ}}}=\underset{i=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{j=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{k=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\frac{V_n(x_n+i,y_n+j,z_n+k)}{{(2\·\mathrm{\δ}+1)}^3},$ 其中i，j，k是临时循环变量；

C03：计算每一个均匀性矩阵V_{Homo_n}的阈值Threshold_{Homo_n}，1≤n≤N，n是当前缩小级数，N是缩小级数的最大值，步骤包括：

C0301生成第n级缩小三维图像V_n的均匀性矩阵V_{Homo_n}的直方图

h_n(r_kn)＝p_kn，其中r_kn是均匀性强度值，p_kn是矩阵中均匀性强度为r_kn的元素个数， 是均匀性矩阵V_{Homo_n}中的最大值；

C0302计算直方图上由0到r_kn区域的高斯拟合函数和其中E(r_kn)是直方图上0到r_kn区域均匀性的均值，Var(r_kn)是直方图上0到r_kn区域均匀性的方差，i是临时循环变量；

C0303计算直方图上由0到r_kn区域的高斯拟合函数和直方图的匹配程度， ${\mathrm{GaussFit}}_n\left(r_{\mathrm{kn}}\right)=\underset{i=0}{\overset{r_{\mathrm{kn}}}{\mathrm{\Σ}}}{h}_n\left(i\right)\·\frac{e^{\frac{{-(i-E\left(r_{\mathrm{kn}}\right))}^2}{2\·\mathrm{Var}\left(r_{\mathrm{kn}}\right)}}}{\sqrt{2\·\mathrm{\πVar}\left(r_{\mathrm{kn}}\right)}},$ i是临时循环变量；

C0304在区域内重复步骤C0302-C0303计算每一个r_kn值的高斯拟合函数和直方图的匹配程度GaussFit_n(r_kn)，将匹配程度的最大值所对应的高斯模型记为E_{Max_n}和Var_{Max_n}，计算阈值 ${\mathrm{Threshould}}_n=E_{\mathrm{Max}\_n}+4\·\sqrt{{\mathrm{Var}}_{\mathrm{Max}\_n}}$

C04：对于投影平面上的任意点的坐标记为^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)其中D₁是X光源到UVW坐标原点ISO的距离；D₂是UVW坐标原点ISO到接收器的距离；u_D是投影面上的点在U坐标轴的坐标值；v_D是投影面上的点在V坐标轴的坐标值；u是虚拟投影面上的对应点在U坐标轴的坐标值；v是虚拟投影面上的对应点在V坐标轴的坐标值；

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

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

C0601设定二维投影上点^UVWP_D(u_D，v_D，-D₂)的灰度值I(u_D，v_D，-D₂)＝0，初始点在W坐标轴的坐标其中s是三维图像对角线的长度 $s=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ I(u_D，v_D，-D₂)是沿着当前射线的总吸收量，u_D表示该直线与投影面交点在U轴上的坐标；v_D表示该直线与虚拟投影面交点在V轴上的坐标，D₂是UVW坐标原点ISO到接收器的距离；

C0602计算当前w坐标对应的点在XYZ坐标系下的坐标取整后的值^XYZP_{l_Z}

C0603设n＝N，n为当前的缩小级数，N为缩小级数的最大值；

C0604将坐标值代入第n级缩小三维图像V_n的均匀性矩阵V_{Homo_n}，并且和对应的Threshold_n进行比较，当 $V_{\mathrm{Homo}\_n}\left(\frac{P_{l\_Z}^{\mathrm{XYZ}}}{2^n}\right)\≤{\mathrm{Threshould}}_n$ 则跳转到步骤C0606，

C0605令n＝n-1，如果n＞0则跳转至步骤C0604；

C0606令I(u_D，v_D，-D₂)＝I(u_D，v_D，-D₂)+2ⁿ·δ·V(^XYZP_{l_Z})且w＝w+2ⁿ·δ，其中I(u_D v_D，-D₂)是沿着当前射线的总吸收量，u_D表示该直线与投影面交点在U轴上的坐标；v_D表示该直线与虚拟投影面交点在V轴上的坐标，D₂是UVW坐标原点ISO到接收器的距离，^XYZP_{l_Z}是直线上点在XYZ坐标系下取整后的坐标；

C0607如果则跳转到步骤C0602，其中s是三维图像对角线的

长度 $s=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2};$

C07：根据Beer’s定理投影平面上对应点的灰度值 $H(u_D,v_D,-D_2)=H_{\max }{e}^{-I(u_D,v_D,-D_2)};$ 其中H_max是投影图像的最大灰度值；

I(u_D，v_D，-D₂)是沿着当前射线的总吸收量，e是自然对数的底数；

C08：对接收器上每一个像素重复步骤C04-C07得到仿真投影DRR。

在上述的论述中^XYZP均代表点在XYZ坐标系的投影，^UVWP均代表点在UVW坐标系的投影，其后括号内的字母均为投影的坐标值。

本发明可以获得如下有益效果：本发明提出的一种基于图像内容均匀性的Raycasting改进算法，在对数据进行预处理时生成用来保存三维数据内容均匀性程度的数据结构，在实际计算光路的过程中根据图像内容的均匀性实现采样点的跳跃，大大提高了DRR生成的速度。

附图说明

图1投影模型以及参数。

图2DRR生成流程图

图3沿X射线积分方法

具体实施方式

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

首先，将三维图像V由DICOM服务器调入工作PC，三维图像V的大小是Width×Height×Depth，其中Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，三维图像V中像素的坐标由(x，y，z)表示。

如图1所示，投影模型由三维图像和投影结构两部分组成，他们分别在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坐标原点到接收器的距离。

如图二所示，是仿真投影生成流程图。生成过程包括如下步骤：第一步，生成三维图像V的缩小图像V₁...V_n...V_N，其中V_n代表第n级的缩小三维图像，其大小为1≤n≤N，n是当前缩小级数，N是缩小级数的最大值，对于大小为512x512x400的三维图像N的最优值为2。第n级缩小三维图像V_n中像素的坐标由(x_n，y_n，z_n)表示，第n级的缩小三维图像V_n中像素和三维图像V中的像素存在对应关系V_n(x_n.y_n.z_n)＝V(2ⁿ·x_n，2ⁿ·y_n，2ⁿ·z_n)。第二步，计算每一级缩小三维图像的均匀性矩阵。以第n级为例，缩小三维图像V_n的均匀性矩阵为V_{Homo_n}，其大小为1≤n≤N，N是缩小级数的最大值，第n级缩小三维图像V_n中坐标为(x_n，y_n，z_n)的像素的均匀性值为 $V_{\mathrm{Homo}\_n}(x_n,y_n,z_n)=\underset{i=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{j=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{k=-\δ}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\frac{{(V_n(x_n+i,y_n+j,z_n+k)-\stackrel{\‾}{V_{x_n{y}_n{z}_n\mathrm{\δ}}})}^2}{{(2\·\mathrm{\δ}+1)}^3},$ 其中δ是半窗长，是第n级缩小三维图像V_n中以像素(x_n，y_n，z_n)为中心，δ为半径覆盖区域像素的平均值 $\stackrel{\‾}{V_{x_n{y}_n{z}_n\mathrm{\δ}}}=\underset{i=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{j=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{k=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\frac{V_n(x_n+i,y_n+j,z_n+k)}{{(2\·\mathrm{\δ}+1)}^3},$ 其中i，j，k是临时循环变量。第三步，计算每一级的阈值。以第n级为例，计算第n级缩小三维图像V_n的均匀性矩阵V_{Homo_n}的阈值Threshold_{Homo_n}，1≤n≤N，n是当前缩小级数，N是缩小级数的最大值，步骤包括：

1生成第n级缩小三维图像V_n的均匀性矩阵V_{Homo_n}的直方图

h_n(r_kn)＝p_kn，其中r_kn是均匀性强度值，p_kn是矩阵中均匀性强度为r_kn的元素个数， 是均匀性矩阵V_{Homo_n}中的最大值；

2计算直方图上由0到r_kn区域的高斯拟合函数和其中E(r_kn)是直方图上0到r_kn区域均匀性的均值，Var(r_kn)是直方图上0到r_kn区域均匀性的方差，i是临时循环变量；

3计算直方图上由0到r_kn区域的高斯拟合函数和直方图的匹配程度， ${\mathrm{GaussFit}}_n\left(r_{\mathrm{kn}}\right)=\underset{i=0}{\overset{r_{\mathrm{kn}}}{\mathrm{\Σ}}}{h}_n\left(i\right)\·\frac{e^{\frac{{-(i-E\left(r_{\mathrm{kn}}\right))}^2}{2\·\mathrm{Var}\left(r_{\mathrm{kn}}\right)}}}{\sqrt{2\·\mathrm{\πVar}\left(r_{\mathrm{kn}}\right)}},$ i是临时循环变量；

4在区域内重复步骤2-3计算每一个r_kn值的高斯拟合函数和直方图的匹配程度GaussFit_n(r_kn)，将匹配程度的最大值所对应的高斯模型记为E_{Max_n}和Var_{Max_n}，计算阈值 ${\mathrm{Threshould}}_n=E_{\mathrm{Max}\_n}+4\·\sqrt{{\mathrm{Var}}_{\mathrm{Max}\_n}}$

如图2所示，二维图像的大小为Width_2DxHeight_2D。为了遍历每一个像素引入临时循环变量i，j。接下来对每一个像素按照图3所示的方法求灰度值的积分。如图3所示，投影平面上点^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点的灰度值是V(X，Y，Z)。则取整后的点^XYZP_{l_Z}的灰度值为V(^XYZP_{l_Z})。对范围 $w\∈\left[\begin{array}{cc}-\frac{s}{2}& \frac{s}{2}\end{array}\right]$ 内射线上所有强度值得积分可以由如下步骤完成：

1设定二维投影上点^UVWP_D(u_D，v_D，-D₂)的灰度值I(u_D，v_D，-D₂)＝0，初始点在W坐标轴的坐标其中s是三维图像对角线的长度 $s=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ I(u_D，v_D，-D₂)是沿着当前射线的总吸收量，u_D表示该直线与投影面交点在U轴上的坐标；v_D表示该直线与虚拟投影面交点在V轴上的坐标，D₂是UVW坐标原点ISO到接收器的距离；

2计算当前w坐标对应的点在XYZ坐标系下的坐标取整后的值^XYZP_{l_Z}

3设n＝N，n为当前的缩小级数，N为缩小级数的最大值；

4将坐标值代入第n级缩小三维图像V_n的均匀性矩阵V_{Homo_n}，并且和对应的Threshold_n进行比较，当 $V_{\mathrm{Homo}\_n}\left(\frac{P_{l\_Z}^{\mathrm{XYZ}}}{2^n}\right)\≤{\mathrm{Threshould}}_n$ 则跳转到步骤6，

5令n＝n-1，如果n＞0则跳转至步骤4；

6令I(u_D，v_D，-D₂)＝I(u_D，v_D ，-D₂)+2ⁿ·δ·V(^XYZP_{l_Z})且w＝w+2ⁿ·δ，其中I(u_D，v_D，-D₂)是沿着当前射线的总吸收量，u_D表示该直线与投影面交点在U轴上的坐标；v_D表示该直线与虚拟投影面交点在V轴上的坐标，D₂是UVW坐标原点ISO到接收器的距离，^XYZP_{l_Z}是直线上点在XYZ坐标系下取整后的坐标；

7如果则跳转到步骤2，其中s是三维图像对角线的长度 $s=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2}$

最后，用Beer’s理论求出该像素的灰度值 $H(u_D,v_D,-D_2)=H_{\max }{e}^{-I(u_D,v_D,-D_2)},$ 其中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 (1)

1.一种基于RayCasting改进算法的仿真投影DRR生成方法，其特征在于，该方法包括如下步骤：

A将三维图像V由DICOM服务器调入工作PC，三维图像V的大小是Width×Height×Depth,其中Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽；

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

B01：以三维图像的中心为坐标原点建立笛卡尔坐标系即三维图像坐标系XYZ，使得X，Y，Z轴分别与三维图像中对应的外平面正交，记坐标原点为O；三维图像V中像素的坐标由(x,y,z)表示，x为像素在X轴上的投影，y为像素在Y轴上的投影,z为像素在Z轴上的投影；

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

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

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

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

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

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

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

C01：计算三维图像V的缩小图像V₁...V_n...V_N，其中V_n代表第n级的缩小三维图像，其大小为1≤n≤N，n是当前缩小级数，N是缩小级数的最大值，第n级缩小三维图像V_n中像素的坐标由(x_n,y_n,z_n)表示，第n级的缩小三维图像V_n中像素和三维图像V中的像素存在对应关系V_n(x_n,y_n,z_n)＝V(2ⁿ·x_n,2ⁿ·y_n,2ⁿ·z_n)；

C02：生成每一个缩小三维图像V_n的均匀性矩阵V_{Homo_n}，其大小为1≤n≤N，N是缩小级数的最大值，第n级缩小三维图像V_n中坐标为(x_n,y_n,z_n)的像素的均匀性值为 $V_{\mathrm{Homo}\_n}(x_n,y_n,z_n)=\underset{j=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{j=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{k=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\frac{{(V_n(x_n+i,y_n+j,z_n+k)-\stackrel{\‾}{V_{x_n,y_n,z_n\mathrm{\δ}}})}^2}{{(2\·\mathrm{\δ}+1)}^3},$ 其中δ是半窗长，是第n级缩小三维图像V_n中以像素(x_n,y_n,z_n)为中心，δ为半径覆盖区域像素的平均值 $\stackrel{\‾}{V_{x_n,y_n,z_n\mathrm{\δ}}}=\underset{i=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{j=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\underset{k=-\mathrm{\δ}}{\overset{\mathrm{\δ}}{\mathrm{\Σ}}}\frac{V_n(x_n+i,y_n+j,z_n+k)}{{(2\·\mathrm{\δ}+1)}^3},$ 其中i,j,k是临时循环变量；C03：计算每一个均匀性矩阵V_{Homo_n}的阈值Threshold_{Homo_n}，1≤n≤N，n是当前缩小级数，N是缩小级数的最大值，步骤包括：

C0301生成第n级缩小三维图像V_n的均匀性矩阵V_{Homo_n}的直方图h_n(r_kn)＝p_kn，其中r_kn是均匀性强度值，p_kn是矩阵中均匀性强度为r_kn的元素个数，0≤r_kn≤是均匀性矩阵V_{Homo_n}中的最大值；

C0302计算直方图上由0到r_kn区域的高斯拟合函数和其中E(r_kn)是直方图上0到r_kn区域均匀性的均值，Var(r_kn)是直方图上0到r_kn区域均匀性的方差，i是临时循环变量；

C0303计算直方图上由0到r_kn区域的高斯拟合函数和直方图的匹配程度， ${\mathrm{GaussFit}}_n\left(r_{\mathrm{kn}}\right)=\underset{i=0}{\overset{r_{\mathrm{kn}}}{\mathrm{\Σ}}}{h}_n\left(i\right)\·\frac{e^{\frac{-{(i-E\left(r_{\mathrm{kn}}\right))}^2}{2\·\mathrm{Var}\left(r_{\mathrm{kn}}\right)}}}{\sqrt{2\·\mathrm{\π}\·\mathrm{Var}\left(r_{\mathrm{kn}}\right)}},$ i是临时循环变量；

C0304在区域0≤r_kn≤内重复步骤C0302-C0303计算每一个r_kn值的高斯拟合函数和直方图的匹配程度GaussFit_n(r_kn)，将匹配程度的最大值所对应的高斯模型记为E_{Max_n}和Var_{Max_n}，计算阈值 ${\mathrm{Threshold}}_n=E_{\mathrm{Max}\_n}+4\·\sqrt{{\mathrm{Var}}_{\mathrm{Max}\_n}}$

C04：对于投影平面上的任意点的坐标记为^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坐标轴的坐标值；

C05：对由坐标为^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的距离；

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

C0601设定二维投影上点^UVWP_D(u_D,v_D,-D₂)的灰度值I(u_D,v_D,-D₂)＝0,初始点在W坐标轴的坐标其中s是三维图像对角线的长度I(u_D,v_D,-D₂)是沿着当前射线的总吸收量，u_D表示该直线与投影面交点在U轴上的坐标；v_D表示该直线与虚拟投影面交点在V轴上的坐标，D₂是UVW坐标原点ISO到接收器的距离；

C0602计算当前w坐标对应的点在XYZ坐标系下的坐标取整后的值^XYZP_{l_Z}

C0603设n＝N，n为当前的缩小级数，N为缩小级数的最大值；

C0604将坐标值代入第n级缩小三维图像V_n的均匀性矩阵V_{Homo_n}，并且和对应的Threshold_n进行比较，当≤Threshold_n则跳转到步骤C0606，

C0605令n＝n-1，如果n>0则跳转至步骤C0604；

C0606令I(u_D,v_D,-D₂)＝I(u_D,v_D,-D₂)+2ⁿ·δ·V(^XYZP_{l_Z})且w＝w+2ⁿ·δ，其中I(u_D,v_D,-D₂)是沿着当前射线的总吸收量，u_D表示该直线与投影面交点在U轴上的坐标；v_D表示该直线与虚拟投影面交点在V轴上的坐标，D₂是UVW坐标原点ISO到接收器的距离，^XYZP_{l_Z}是直线上点在XYZ坐标系下取整后的坐标；

C0607如果则跳转到步骤C0602，其中s是三维图像对角线的

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

C08：对接收器上每一个像素重复步骤C04-C07得到仿真投影DRR。