[发明专利]一种基于RayCasting改进算法的仿真投影DRR生成方法

摘要

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

主权项

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；其公式表示如下：TUVWXYZ=cosβ0sinβ00100-sinβ0cosβ0000110000cosα-sinα00sinαcosα00001cosφ-sinφ00sinφcosφ0000100001100Tu010Tv00100001]]>α为投影模型坐标系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)的像素的均匀性值为VHomo_n(xn,yn,zn)=Σi=-δδΣj=-δδΣk=-δδ(Vn(xn+i,yn+j,zn+k)-Vxnynznδ‾)2(2·δ+1)3,]]>其中δ是半窗长，是第n级缩小三维图像V_n中以像素(x_n，y_n，z_n)为中心，δ为半径覆盖区域像素的平均值Vxnynznδ‾=Σi=-δδΣj=-δδΣk=-δδVn(xn+i,yn+j,zn+k)(2·δ+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区域的高斯拟合函数和直方图的匹配程度，GaussFitn(rkn)=Σi=0rknhn(i)·e-(i-E(rkn))22·Var(rkn)2·πVar(rkn),]]>i是临时循环变量；C0304在区域内重复步骤C0302-C0303计算每一个r_kn值的高斯拟合函数和直方图的匹配程度GaussFit_n(r_kn)，将匹配程度的最大值所对应的高斯模型记为E_{Max_n}和Var_{Max_n}，计算阈值Threshouldn=EMax_n+4·VarMax_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)其中u=uDD1D1+D2v=vDD1D1+D2;]]>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；其中XYZTUVW是坐标系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=Width2+Height2+Depth2,]]>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进行比较，当VHomo_n(Pl_ZXYZ2n)≤Threshouldn]]>则跳转到步骤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=Width2+Height2+Depth2;]]>C07：根据Beer’s定理投影平面上对应点的灰度值H(uD,vD,-D2)=Hmaxe-I(uD,vD,-D2);]]>其中H_max是投影图像的最大灰度值；I(u_D，v_D，-D₂)是沿着当前射线的总吸收量，e是自然对数的底数；C08：对接收器上每一个像素重复步骤C04-C07得到仿真投影DRR。