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

賈曉未，魏 嵬，賈克斌

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

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

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

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

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

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

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

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

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

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

1.2 傳統(tǒng)方法

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

1.3 投影模型

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

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

通過

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

2 基于圖元的加速方法

2.1 方法描述

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

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

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

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

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

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

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

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

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

整理后得

由克萊默法則，有

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

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

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

3 基于GPU的并行加速

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

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

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

當求交過程結(jié)束后，需要進行積分的計算，以每個像素為線程，過程如下:

(1)對交點位置進行排序。

(2)依次判斷，當兩個相鄰的交點分別為入射和出射的狀態(tài)時，計算兩點間積分，累加到該像素的積分值。

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

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

改進算法完整流程如圖5所示。

圖5 改進的MeshDRR算法流程圖Fig.5 Improved MeshD RRGeneration Method

4 實驗結(jié)果與比較

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

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

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

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

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

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

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

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

4 結(jié)束語

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

［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.