[發明專利]一種基于RayCasting改進算法的仿真投影DRR生成方法

有效

• 申請號: 201210109477.4
• 申請日: 2012-04-13
• 公開（公告）號: CN102663803A
• 公開（公告）日: 2012-09-12
• 發明（設計）人: 賈克斌;魏嵬;賈曉未
• 申請（專利權）人: 北京工業大學
• 主分類號: G06T15/06
• 分類號: G06T15/06;G06T15/60
• 代理公司: 北京思海天達知識產權代理有限公司 11203
• 代理人: 魏聿珠

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。