一種基于GPU的圖元網(wǎng)狀結(jié)構(gòu)DRR并行加速算法

賈曉未，魏 嵬，賈克斌

(1.CSEDepartment，Universityat Buf falo，SUNY14260;2.北京工業(yè)大學(xué) 電子信息與控制學(xué)院，北京 100124)

計(jì)算機(jī)輔助導(dǎo)航手術(shù)(CAS－Computer Assisted Surgery)系統(tǒng)是計(jì)算機(jī)圖像處理及三維可視化技術(shù)、機(jī)器人技術(shù)、空間三維定位導(dǎo)航技術(shù)和臨床手術(shù)等多學(xué)科技術(shù)結(jié)合的產(chǎn)物。其目的主要是借用計(jì)算機(jī)技術(shù)及設(shè)備跟蹤模擬手術(shù)中涉及的每個(gè)過(guò)程。計(jì)算機(jī)輔助手術(shù)導(dǎo)航系統(tǒng)在骨科創(chuàng)傷手術(shù)微創(chuàng)化進(jìn)程中有著極其重要的意義。計(jì)算機(jī)輔助骨外科導(dǎo)航手術(shù)(CAOS－Computer Assisted Orthopedic Surgery)通過(guò)計(jì)算機(jī)輔助的術(shù)前規(guī)劃，模擬手術(shù)，術(shù)中追蹤手術(shù)器械，引導(dǎo)手術(shù)以及術(shù)后評(píng)估，使得骨科創(chuàng)傷手術(shù)更安全，更快速，更準(zhǔn)確。其中二維三維配準(zhǔn)是系統(tǒng)開發(fā)中的一項(xiàng)關(guān)鍵技術(shù)。目前二維三維圖像配準(zhǔn)還處于研究階段［1］，還沒有穩(wěn)定的方法可以應(yīng)用于臨床。在現(xiàn)今存在的配準(zhǔn)技術(shù)中，基于DRR(DigitallyReconstructed Radiograph)的配準(zhǔn)方法使用廣泛。其中DRR生成也是配準(zhǔn)過(guò)程中運(yùn)算時(shí)間最長(zhǎng)的一個(gè)單元，作為配準(zhǔn)過(guò)程的中間介質(zhì)，DRR的生成速度和質(zhì)量往往對(duì)配準(zhǔn)的效率和精度有至關(guān)重要的影響。

目前在實(shí)際二維三維配準(zhǔn)應(yīng)用中，除了CT圖與2DX 光片直接利用圖片信息的配準(zhǔn)［2－4］，還有很多應(yīng)用通過(guò)植入器件的在3D和2D空間中的配準(zhǔn)實(shí)現(xiàn)不同維度下的位置匹配［5－6］，植入器件的方法由于其較為明顯的外部特征在目前的手術(shù)應(yīng)用中也十分常見。不同于CT數(shù)據(jù)的Volume結(jié)構(gòu)，植入器件往往以網(wǎng)狀結(jié)構(gòu)(Mesh)存儲(chǔ)，這種結(jié)構(gòu)相比于Volume結(jié)構(gòu)具有空間的不規(guī)則性，因此也給DRR的生成帶來(lái)了很大的難度，傳統(tǒng)的方法在遍歷時(shí)需要消耗大量時(shí)間。本文在傳統(tǒng)的RayCasting 算法［7－8］基礎(chǔ)上，提出了一種基于圖元的加速方法，并且給出了針對(duì)空間平面與直線的求交問(wèn)題和光線篩選問(wèn)題的優(yōu)化方法。實(shí)驗(yàn)結(jié)果表明，這一方法大大降低了網(wǎng)狀結(jié)構(gòu)DRR生成的時(shí)間消耗，同時(shí)該方法也具有可并行性。

1 網(wǎng)狀結(jié)構(gòu)及傳統(tǒng)生成方法

1.1 網(wǎng)狀(Mesh)結(jié)構(gòu)

網(wǎng)狀結(jié)構(gòu)的數(shù)據(jù)由眾多圖元組成，旨在對(duì)圖元的空間結(jié)構(gòu)進(jìn)行存儲(chǔ)。對(duì)于植入器件，圖元通常呈三角形狀，通過(guò)記錄三個(gè)頂點(diǎn)在三維空間坐標(biāo)系下的坐標(biāo)保存每個(gè)圖元的空間位置。完整的植入模型數(shù)據(jù)通常由大量的圖元組成。

如圖1所示，網(wǎng)狀結(jié)構(gòu)中每個(gè)三角形片層以其頂點(diǎn)A，B，C的空間坐標(biāo)依次存儲(chǔ)，向量AB×AC表示了該平面的法線方向，指向物體外側(cè)。網(wǎng)狀結(jié)構(gòu)模型由這樣的片層結(jié)構(gòu)組成其表面，通常其內(nèi)部認(rèn)為是勻質(zhì)結(jié)構(gòu)，因此DRR中像素值與對(duì)應(yīng)光線的射入光路距離直接相關(guān)。

圖1 網(wǎng)狀結(jié)構(gòu)圖元示意圖Fig.1 Triang leprimitives tructure

圖2為一個(gè)植入器件的網(wǎng)狀結(jié)構(gòu)模型。對(duì)一個(gè)勻質(zhì)結(jié)構(gòu)，通常在保證失真程度較小的前提下將其空間結(jié)構(gòu)以網(wǎng)狀形式存儲(chǔ)，用網(wǎng)狀結(jié)構(gòu)中每個(gè)三角形片元頂點(diǎn)的空間坐標(biāo)記錄其空間位置，以此保證整個(gè)模型的外表面結(jié)構(gòu)被較為完整地保存。

圖2 投影模型示意圖Fig.2 Projection Model

1.2 傳統(tǒng)方法

傳統(tǒng)的Ray Casting算法通過(guò)遍歷每個(gè)像素，并且計(jì)算每條由像素投射出的光線被遮擋部分的密度積分，最終合成整幅DRR圖像。對(duì)于網(wǎng)狀結(jié)構(gòu)，與Volume結(jié)構(gòu)不同的是，由于不能根據(jù)光線的斜率直接計(jì)算出相關(guān)的交點(diǎn)，往往需要遍歷所有的三角形面，并依次求取交點(diǎn)然后計(jì)算積分值，最終合成DRR圖像。通過(guò)搭建如八叉樹等數(shù)據(jù)結(jié)構(gòu)，可以提升遍歷的速度。但是由于網(wǎng)狀結(jié)構(gòu)對(duì)于投影模型不確定的空間特性，搭建這樣的數(shù)據(jù)結(jié)構(gòu)具有一定的困難。在傳統(tǒng)RayCasting算法上，本文將提出一種基于圖元的加速算法，很大程度上減少了計(jì)算量，且算法具有可并行性。

1.3 投影模型

本文引入一種投影模型，如圖3所示。其中UVW為投影模型坐標(biāo)系，XYZ為數(shù)據(jù)模型坐標(biāo)系。Ru，Rv，Rw 表示繞 U，V，W 軸的旋轉(zhuǎn)角，Tu，Tv表示在UV平面內(nèi)的平移參數(shù)，D1，D2分別表示光源距ISO中心和ISO中心距接收屏的距離。對(duì)應(yīng)的轉(zhuǎn)換矩陣為:

圖3 投影模型示意圖Fig.3 Projection Model

通過(guò)

可將UVW坐標(biāo)系下坐標(biāo)轉(zhuǎn)至XYZ坐標(biāo)系。

2 基于圖元的加速方法

2.1 方法描述

本文以每個(gè)圖元，即三角形片作為一個(gè)單元，對(duì)于由像素投射出的光線進(jìn)行遍歷，計(jì)算光線與三角形平面的交點(diǎn)并記錄在以光線序號(hào)為下標(biāo)的結(jié)構(gòu)中，完成所有圖元的計(jì)算后通過(guò)每條光線的積分值合成DRR圖像。

實(shí)際的植入器件往往由大量圖元片層構(gòu)成，因此基于圖元的算法具有更強(qiáng)的并行性。在另一方面，對(duì)于某一三角面圖元，可以通過(guò)光源與接收屏的位置進(jìn)行光線篩選，縮小需要遍歷的光線范圍，進(jìn)而很大程度上減少了每個(gè)圖元的計(jì)算任務(wù)量，算法過(guò)程如圖3所示。

圖3 基于圖元的Ray Casting方法Fig.3 Primitive-based Ray Casting Method

由于真實(shí)數(shù)據(jù)中圖元數(shù)量較大，因此圖3中第一個(gè)對(duì)于圖元的遍歷可以通過(guò)硬件并行計(jì)算，能夠顯著減少時(shí)間消耗。對(duì)于光線范圍的篩選，首先將三角面三個(gè)頂點(diǎn)坐標(biāo)進(jìn)行式(2)的逆變換，由XYZ坐標(biāo)系轉(zhuǎn)換至UVW坐標(biāo)系，然后在UVW坐標(biāo)系下計(jì)算由光源，即(O，O，Dl)點(diǎn)發(fā)出的經(jīng)過(guò)這三個(gè)頂點(diǎn)的三條直線與平面W=－D2，即接收屏的交點(diǎn)，分別取這三個(gè)交點(diǎn)在U方向和V方向上的最小值和最大值，以此得到所需遍歷的光線的最小矩形范圍區(qū)域。

2.2 空間中求交點(diǎn)的優(yōu)化方法

在網(wǎng)狀結(jié)構(gòu)生成DRR的過(guò)程中，核心問(wèn)題在于判斷空間中直線是否與三角面相交，并且計(jì)算交點(diǎn)。處理這個(gè)問(wèn)題消耗的時(shí)間將直接影響整個(gè)算法的效率。

實(shí)際應(yīng)用中往往通過(guò)先求直線與三角形所在平面的交點(diǎn)，再判斷該交點(diǎn)是否位于三角形平面內(nèi)。由于判斷同一平面中的一點(diǎn)是否位于三角形內(nèi)相對(duì)容易，因此這種方法也被廣為接受。然而，由于網(wǎng)狀結(jié)構(gòu)中存在大量十分微小的三角面，而且光線數(shù)目眾多，因此這種方法往往在計(jì)算交點(diǎn)上花費(fèi)了過(guò)多的時(shí)間。

本文采用類似文獻(xiàn)［9］的方法求解空間三角面與直線的交點(diǎn)。首先以

分別表示空間三角面與直線，式中:(u，v)為一點(diǎn)在三角形內(nèi)重心坐標(biāo)系下的坐標(biāo);O為光源點(diǎn);D為步長(zhǎng)。聯(lián)立方程(3)和方程(4)有

整理后得

由克萊默法則，有

式中:E1=B －A，E2=C －A，T=O －A，P=(D ×E2)，Q=(T×E1)。交點(diǎn)在三角形內(nèi)的充要條件為 u≥0，v≥0，u+v≤1。

為了對(duì)光線的積分值進(jìn)行計(jì)算，需要判斷光線對(duì)于該三角面為入射狀態(tài)還是出射狀態(tài)。在1.1節(jié)中提到，E1×E2表示法線方向，并指向物體外側(cè)，因此通過(guò)判斷(L=D·(E1×E2))是否大于0來(lái)確定出射點(diǎn)和入射點(diǎn)，可以得到入射或出射的信息。

在圖3中最后對(duì)于光線的循環(huán)中，當(dāng)按序排列的交點(diǎn)列中相鄰的兩點(diǎn)依次為入射和出射時(shí)，將兩點(diǎn)間的積分值累加到整條光線的積分值，最終合成整幅DRR圖像。此方法計(jì)算過(guò)程中考慮了一些多次計(jì)算的中間變量，同時(shí)可以根據(jù)變量的值判斷三角片元與光線是否相交，若不相交則不用準(zhǔn)確計(jì)算交點(diǎn)位置，因此相對(duì)于傳統(tǒng)方法節(jié)省了不必要的計(jì)算時(shí)間。

3 基于GPU的并行加速

基于圖元的Mesh DRR生成算法具有較好的并行性，且對(duì)于光線的篩選優(yōu)化可以有效減少每個(gè)線程的計(jì)算任務(wù)量。

在CUDA架構(gòu)下，顯示芯片執(zhí)行時(shí)的最小運(yùn)算單元是工作線程(Thread)。多個(gè)工作線程組成一個(gè)塊(Block)。多個(gè)塊又組成一個(gè)柵格(Grid)。每一個(gè)柵格對(duì)應(yīng)著一個(gè)設(shè)備程序即Kernel。對(duì)于本文中提出的基于圖元的網(wǎng)狀結(jié)構(gòu)模型DRR生成算法，以針對(duì)每個(gè)圖元的計(jì)算量為GPU中的一個(gè)線程并行計(jì)算。

由于CUDA對(duì)資源訪問(wèn)的限制，不同塊之間的內(nèi)存是不能共享的，不同塊之間的通信也受到一定限制。由于不同線程對(duì)同一像素值進(jìn)行改寫時(shí)會(huì)發(fā)生寫沖突，本文中主要采用重啟Kernel函數(shù)的方法實(shí)現(xiàn)線程同步，同時(shí)，設(shè)定一定容量的暫存空間保存中間結(jié)果，以減少重啟Kernel函數(shù)的次數(shù)，進(jìn)而減少了因重啟Kernel函數(shù)帶來(lái)的開銷。首先建立以圖元編號(hào)為下標(biāo)的結(jié)構(gòu)，在Kernel函數(shù)中每個(gè)線程將交點(diǎn)信息及光線編號(hào)存入對(duì)應(yīng)的位置中，當(dāng)顯存占滿時(shí)，退出Kernel函數(shù)，將數(shù)據(jù)導(dǎo)入以光線編號(hào)為下標(biāo)的結(jié)構(gòu)中，然后再啟動(dòng)Kernel函數(shù)，直至光線遍歷結(jié)束。這種方法的缺陷在于，兩種儲(chǔ)存結(jié)構(gòu)的轉(zhuǎn)換需要消耗大量的時(shí)間。因此，本文中將這種轉(zhuǎn)換以Kernel函數(shù)的形式在GPU中執(zhí)行，以對(duì)目標(biāo)結(jié)構(gòu)每個(gè)光線編號(hào)對(duì)應(yīng)的位置進(jìn)行存儲(chǔ)的任務(wù)作為一個(gè)線程，以此減少時(shí)間消耗。

當(dāng)求交過(guò)程結(jié)束后，需要進(jìn)行積分的計(jì)算，以每個(gè)像素為線程，過(guò)程如下:

(1)對(duì)交點(diǎn)位置進(jìn)行排序。

(2)依次判斷，當(dāng)兩個(gè)相鄰的交點(diǎn)分別為入射和出射的狀態(tài)時(shí)，計(jì)算兩點(diǎn)間積分，累加到該像素的積分值。

(3)若不滿足入射出射的狀態(tài)，則說(shuō)明該交點(diǎn)位于圖像的棱角處，不需要進(jìn)行累加積分，尋找下一組滿足要求的交點(diǎn)。

4)每個(gè)線程判斷完所有的交點(diǎn)得到最終每個(gè)像素的積分值，合成整幅DRR圖像。

改進(jìn)算法完整流程如圖5所示。

圖5 改進(jìn)的MeshDRR算法流程圖Fig.5 Improved MeshD RRGeneration Method

4 實(shí)驗(yàn)結(jié)果與比較

所用實(shí)驗(yàn)數(shù)據(jù)為圖6所示的植入模型的真實(shí)樣例。

圖6 網(wǎng)狀結(jié)構(gòu)植入模型示意圖Fig.6 Implant Model

實(shí)驗(yàn)環(huán)境為IntelXeonCPUE5405@2.00 GHz2.00GHz，3.25GBofRAM，顯卡 Quadro FX570，顯存 256MB。用于實(shí)驗(yàn)的程序在WindowXP平臺(tái)下開發(fā)，由 VS2008編寫，其中GPU加速部分基于Cuda1.0。數(shù)據(jù)為35378片的植入模型壓縮數(shù)據(jù)，實(shí)驗(yàn)生成DRR圖像，圖像分辨率為200×250。

以本文算法生成DRR圖像如圖7所示。

圖7 植入模型不同角度DRR生成圖Fig.7 DRR of Implant Modelin Dif ferent Views

實(shí)驗(yàn)隨機(jī)產(chǎn)生Ru，Rv，Rw50次，并對(duì)時(shí)間消耗求平均值。表1為傳統(tǒng)算法與改進(jìn)算法的時(shí)間消耗。

表1 傳統(tǒng)算法與改進(jìn)算法的時(shí)間消耗Table1 Time Coston Classic Method and Accelerated Methods

由表1可以看出，經(jīng)過(guò)2.2節(jié)中所提到的交點(diǎn)求取方法優(yōu)化后，時(shí)間消耗有了很好的改善。由于本文采用了基于圖元的方法，因此得以在此基礎(chǔ)上針對(duì)每個(gè)圖元做光線的篩選。由表1可以看出，這一步改進(jìn)對(duì)時(shí)間消耗的改進(jìn)更為顯著，原因主要在于，在真實(shí)數(shù)據(jù)中，存在大量的光線與圖元沒有相交。經(jīng)過(guò)GPU加速的算法由于各個(gè)線程的并行計(jì)算從而使得時(shí)間消耗有明顯減少。

4 結(jié)束語(yǔ)

由實(shí)驗(yàn)數(shù)據(jù)可以看出，這種基于圖元大大減少了時(shí)間消耗，同時(shí)可以保證圖片質(zhì)量。與傳統(tǒng)RayCasting的并行性［10］相比，這種方法同樣適用于并行運(yùn)算，以每個(gè)圖元作為并行計(jì)算的線程，另外光線篩選的方法顯著降低了每個(gè)線程的計(jì)算任務(wù)量。因此這種方法有其應(yīng)用價(jià)值，同時(shí)為基于植入模型的二維三維配準(zhǔn)在性能上的提升奠定了基礎(chǔ)。

［1］Markelj P，Tomazevic D，Likar B，et al.A review of 3D/2D registration methods for image－guided interventions［J］.Medical Image Analysis，2010，16(3):642－661.

［2］Lemieux L，Jagoe R，F(xiàn)ish D R，et al.A patient to computed tomography image registration method based on digitally reconstructed radiograph［J］.Med.Phys，1994，21(11):1749－1760.

［3］Penny G P，Weese J，Little J A，et al.A comparison of similarity measures for use in 2－D－3－D medical image registration［J］.IEEE Transactions on Medical Imaging，1998，17(4):586－595.

［4］Jacquet W，Nyssen E，Bottenberg P，et al.2D image registration using focused mutual information for application in dentistry［J］.Computers in Biology and Medicine，2009，39(6):545－553.

［5］Mohamed R Mahfouz，William A Hoff，Richard D Komistek，et al.A robust method for registration of three－dimensional knee implant models to two－dimensional fluoroscopy images［J］.IEEE Transactions on Medical Imaging，2003，22(12):1561－1574.

［6］Mohamed R Mahfouz，William A Hoff，Richard D Komistek，et al.Effect of segmentation errors on 3D－to－2D registration of implant models in X－ray images［J］.Journal of Biomechanics，2005，38(2):229－239.

［7］James F Blinn.Light reflection functions for simulation of clouds and dusty surfaces［J］.Computer Graphics，1982，16(3):21－29.

［8］Levoy M.Display of surfaces from volume data［D］.Chapel Hill，North Carolina:Department of Computer Science，University of North Carolina at Chapel Hill，1989.

［9］Prosolvia Clarus A B，Ben Trumbore.Fast，minimum storage ray－triangle intersection［J］.Journal of Graphics，GPU＆ Game Tools，1997，2(1):21－28.

［10］Ruijters D，Ter Haar Romeny B M，Suetens P.GPU－accelerated digitally reconstructed radiographs［C］//BioMED(08 Proceedings of the Sixth IASTED International Conference on Biomedical Engineering)，Canada，2008，601:431－435.