專利名稱：一種基于顯卡并行計算的三維醫(yī)學圖像的快速分割方法
技術(shù)領(lǐng)域：
本發(fā)明屬于醫(yī)學圖像處理、醫(yī)學圖像分割技術(shù)領(lǐng)域，特別涉及基于顯卡并行計算的三維醫(yī)學圖像的快速分割方法。

背景技術(shù)：
在可視化醫(yī)療手術(shù)中，人體器官及病灶的三維構(gòu)形技術(shù)為醫(yī)生進行手術(shù)帶來了巨大的便利，而三維構(gòu)形的前提則是從原始的，由大量二維醫(yī)學圖像片層構(gòu)成的三維圖像(例如CT、MRI圖像)中將目標器官或是腫瘤圖像準確的分割出來。這種分割在醫(yī)療圖像處理領(lǐng)域中一直是一個非常困難的問題，尤其是在需要兼顧分割的準確性和處理的快速性的前提下。
過去從醫(yī)療圖像中提取各種器官的圖像都是靠人工在每個圖像片層中進行邊界勾畫來完成的，而該過程的工作量往往需要耗時數(shù)個小時。因此為了減少時間開銷，一些學者開始研究半自動以及自動的器官圖像分割方法。例如韓國的Lee等人在2007年發(fā)表的論文“Efficient liver segmentation using a level-set method with optimal detectionof the initial liver boundary from level-set speed images”中采用了結(jié)合邊界信息的改進的水平集方法來進行肝臟圖像的分割，該方法結(jié)果雖然準確，但需要較多的時間開銷處理一張512x512像素的CT片層需要3.35秒來完成。這在數(shù)據(jù)量較大時是無法滿足臨床的實時性要求的。
Graph-Cuts是另一種可以達到很高分割精度的方法，且這種方法近年來在包括醫(yī)療圖像分割在內(nèi)的很多領(lǐng)域內(nèi)都有著廣泛的應(yīng)用。該方法的提出者Yuri Boykov在2006年發(fā)表的論文“Graph-Cuts and Efficient N-D image segmentation”中即利用其對人體肝臟圖像進行了分割并獲得了不錯的效果。
Graph-Cuts根據(jù)原始輸入的圖像構(gòu)造出一幅對應(yīng)的網(wǎng)絡(luò)流圖，該網(wǎng)絡(luò)流圖由原始圖像像素(二維情況下)或是三維體素(三維情況下)轉(zhuǎn)化的節(jié)點，額外添加的起始端點與終止端點，連接節(jié)點與節(jié)點及節(jié)點與端點間的帶有權(quán)值的邊構(gòu)成；連接節(jié)點與節(jié)點的邊的權(quán)值用于反映輸入圖像中的邊界信息相鄰的兩個像素(或三維體素)的差值越大，則該方法令這兩個像素(或三維體素)轉(zhuǎn)化的節(jié)點間連邊的權(quán)值越小，因為從直觀上講，兩個像素(或三維體素)間的差值大，則這兩個像素(或三維體素)更可能位于分割的邊界上(即被分別分割到目標和背景中)；連接節(jié)點與端點的邊的權(quán)值用于反映輸入圖像中的區(qū)域信息該方法令添加的起始端點和終止端點分別代表分割的背景和目標，節(jié)點的像素(或三維體素)值越接近于預(yù)先獲知的分割目標的取值情況，則該節(jié)點與終止端點間連邊的權(quán)值越小，而節(jié)點的像素(或三維體素)值越接近于預(yù)先獲知的背景的取值情況，則該節(jié)點與起始端點間連邊的權(quán)值越小。Graph-Cuts方法即在上述規(guī)定下計算出網(wǎng)絡(luò)流圖中由起始端點到終止端點的最小割切(網(wǎng)絡(luò)流圖中從起點到終點的割切是指邊的這樣一個集合如果從網(wǎng)絡(luò)流圖中刪去切割中包含的所有邊，則網(wǎng)絡(luò)流圖中將無法在余下的邊中找到一條從起點到終點的路徑。最小割切則是指網(wǎng)絡(luò)流圖的所有割切中包含的各邊權(quán)值之和最小的割切)，依據(jù)上述邊權(quán)值的計算方式，那些差值大的像素(或三維體素)轉(zhuǎn)化的節(jié)點間的邊、像素(或三維體素)值更接近于分割目標的節(jié)點與終止端點間的連邊以及像素(或三維體素)值更接近于背景的節(jié)點與起始端點間的連邊就更可能出現(xiàn)在最小分割中。在Graph-Cuts方法中計算得到的最小割切中，那些連接節(jié)點與節(jié)點間的邊就對應(yīng)著原始圖像中分割的邊界(即每條這樣的邊連接的兩個節(jié)點所對應(yīng)的像素(或三維體素)分別位于分割的目標和背景中)。然而Graph-Cuts算法中盡管給出了網(wǎng)絡(luò)流圖中邊權(quán)值的規(guī)定方式，但沒有對這些權(quán)值的具體計算給出形式化的表述，實際的計算公式需要根據(jù)對具體分割目標的分析和反復(fù)的實驗的總結(jié)才能最終確定。
在圖論中，計算網(wǎng)絡(luò)流圖的最小割切的問題等價于計算該網(wǎng)絡(luò)流圖上的最大流的問題。形象的講，把網(wǎng)絡(luò)流圖的起始端點看作是一個油料的輸送站，終止端點看作接收站，各節(jié)點看作中轉(zhuǎn)站，而所有邊都看作是輸油管道，邊的權(quán)值看作是輸油管道的容量，那么最大流問題就是需要計算出一組輸油路徑，使得盡可能多的油能夠從輸油站到達接收站，而網(wǎng)絡(luò)流圖的最大流的流量(到達接收站的油料數(shù)量)就等于最小割切中各邊的權(quán)值之和，因此最小分割中的邊可以通過在最大流中尋找流量等于權(quán)值(即在輸送過程中該管道被油充滿)的邊來獲得。Push-Relabel方法是一種常用的計算網(wǎng)絡(luò)流圖的方法，該方法初始時允許一些節(jié)點處流入的流量之和大于流出的流量之和，然后對所有這樣的節(jié)點反復(fù)進行push和relabel這兩種基本操作將節(jié)點多出的流量轉(zhuǎn)移到終止端點(接收站)或是返回給起始端點(輸油站)，當所有節(jié)點的流入流量等于流出流量之時，網(wǎng)絡(luò)流圖就達到了最大流的狀態(tài)，也就等價于得到了該網(wǎng)絡(luò)流圖的最小割切。然而該方法的運算復(fù)雜度依然較高，對于n個節(jié)點和m條邊構(gòu)成的網(wǎng)絡(luò)流圖，Push-Relabel方法的時間復(fù)雜度為O(nm2)，在傳統(tǒng)的CPU基礎(chǔ)上實現(xiàn)的最大流算法對于一幅尺寸為170x170x144的三維醫(yī)學圖像需要約30秒來完成計算，這也很難滿足臨床應(yīng)用的實時要求。
近年來，圖形處理單元(GPU)作為一種經(jīng)濟且高效的并行處理設(shè)備而受到了研究界的重視。一塊當前的GPU可以包含多達數(shù)百個的流處理單元，因此其非常適合解決并行的計算問題。Dixit等人在2005年發(fā)表的論文“Gpu-cutsCombinatorial optimization，graphic processing units and adaptive object extraction”中即給出了一種基于早期GPU的Graph-Cuts分割方法的并行實現(xiàn)，然而受限于當時的GPU技術(shù)，該方法的實際運行效率比傳統(tǒng)方法還要慢一些。更新一些的研究則利用NVIDIA公司提出的一種通用并行計算架構(gòu)CUDA來完成Graph-Cuts的并行化處理。在Hussein等人在2007年發(fā)表的“OnImplementing Graph Cuts on CUDA”和Vineet等人在2008年發(fā)表的論文“CUDA CutsFastGraph Cuts on the GPU”中，他們都使用了這一架構(gòu)實現(xiàn)了針對二維圖像的CUDA風格的Graph-Cuts方法，然而他們沒有考慮原始圖像規(guī)模較大時節(jié)省存儲空間的問題，因此并不適用于對存儲空間要求更為嚴格的三維空間下的應(yīng)用。


發(fā)明內(nèi)容
本發(fā)明的目的是為克服已有技術(shù)的不足之處，提出一種基于顯卡并行計算的三維醫(yī)學圖像的快速分割方法，該方法通過同時考量圖像中的區(qū)域信息和邊界信息以獲得最優(yōu)的分割結(jié)果，并利用CUDA技術(shù)進行方法的并行加速，能夠在保持高處理精度的前提下大幅提高分割方法的運行速度。
本發(fā)明提出的一種基于顯卡并行計算的三維醫(yī)學圖像的快速分割方法，其特征在于，該方法采用基于CUDA技術(shù)加速的Graph-Cuts方法對待分割的原始三維醫(yī)學圖像進行分割處理；包括以下步驟 1)將待分割的原始三維醫(yī)學圖像初始化為一個網(wǎng)絡(luò)流圖； 2)在網(wǎng)絡(luò)流圖上采用改進的最大流方法進行分割處理，得到三維醫(yī)學圖像的分割結(jié)果； 上述步驟1)將待分割的原始三維醫(yī)學圖像初始化為一個網(wǎng)絡(luò)流圖，具體包括以下步驟 11)首先根據(jù)原始三維醫(yī)學圖像構(gòu)造出由節(jié)點、端點及有權(quán)邊組成的網(wǎng)絡(luò)流圖的基本架構(gòu)； 12)根據(jù)網(wǎng)絡(luò)流圖的基本構(gòu)架構(gòu)造網(wǎng)絡(luò)流圖的邏輯結(jié)構(gòu)為每一個節(jié)點和端點分配一個整數(shù)值作為該節(jié)點的標號，設(shè)初始時端點S的標號值為N，其它所有節(jié)點和端點的標號值為0；同時為每條有權(quán)邊分配一個絕對值不超過該邊容量值的流量值，初始時設(shè)所有有權(quán)邊的流量值均為0；即構(gòu)成網(wǎng)絡(luò)流圖的邏輯結(jié)構(gòu)； 13)構(gòu)造網(wǎng)絡(luò)流圖的存儲結(jié)構(gòu)從而獲得完成實際存儲的網(wǎng)絡(luò)流圖 所述存儲結(jié)構(gòu)為基于節(jié)點的由兩組五個三維數(shù)組組成，將網(wǎng)絡(luò)流圖中的所有的有權(quán)邊的容量值存儲于第一組的五個三維數(shù)組CoEX，CoEY，CoEZ，CoES和CoET中，將網(wǎng)絡(luò)流圖中的流量值存儲于第二組的五個三維數(shù)組FoEX，F(xiàn)oEY，F(xiàn)oEZ，F(xiàn)oES和FoET中，通過節(jié)點在原始三維醫(yī)學圖像中的坐標來獲得與其相連接的有權(quán)邊的索引；這樣就得到了用于實際存儲的網(wǎng)絡(luò)流圖； 上述步驟2)在實際存儲的網(wǎng)絡(luò)流圖上采用改進的最大流方法進行分割處理，具體包括以下步驟 21)對實際存儲的網(wǎng)絡(luò)流圖初始化活躍節(jié)點作為分割的啟動條件，該活躍節(jié)點為流入該節(jié)點的流量值之和大于流出該節(jié)點的流量值之和的節(jié)點；令網(wǎng)絡(luò)流圖中所有節(jié)點與端點S間的有權(quán)邊的流量值等于該邊的容量值，使網(wǎng)絡(luò)流圖中所有的節(jié)點均成為活躍節(jié)點； 22)對實際存儲的網(wǎng)絡(luò)流圖進行基于CUDA實現(xiàn)的兩階段并行relabel操作 首先僅對偶節(jié)點并行進行relabel操作，偶節(jié)點為節(jié)點三維坐標(a，b，c)滿足a+b+c；然后對余下的奇節(jié)點并行進行relabel操作；以更新所有活躍節(jié)點的標號 23)對完成relabel操作后的網(wǎng)絡(luò)流圖進行基于CUDA實現(xiàn)的兩階段并行push操作首先偶節(jié)點并行進行push操作，然后對余下的奇節(jié)點并行進行push操作；以將活躍節(jié)點的超額流量轉(zhuǎn)移到一個標號值小一的相鄰節(jié)點或端點上 24)當活躍節(jié)點的超額流量值被更新為0，則該活躍節(jié)點就變?yōu)槠胀ü?jié)點，檢查網(wǎng)絡(luò)流圖中是否還存在活躍節(jié)點；如果依然存在，則回轉(zhuǎn)到22)步驟繼續(xù)進行relabel和push操作；否則若網(wǎng)絡(luò)流圖中不再存在活躍節(jié)點，網(wǎng)絡(luò)流圖中將得到從S到T的最大流，最大流對應(yīng)的最小割切中容量值等于流量值的邊即構(gòu)成了分割目標器官與背景的邊界，即得到了分割的結(jié)果。
本發(fā)明的技術(shù)特點及有益效果 本發(fā)明在實現(xiàn)這一分割方法時的創(chuàng)新點是針對數(shù)據(jù)量較大的三維醫(yī)學圖像專門設(shè)計了能夠快速訪問的基于節(jié)點的網(wǎng)絡(luò)流圖存儲結(jié)構(gòu)，并且將Graph-Cuts方法中的Push操作和Relabel操作都進行了CUDA架構(gòu)下的并行化處理以提高運行速度，同時設(shè)計了相鄰節(jié)點間更新數(shù)據(jù)時無需額外存儲消耗的沖突避免的機制。這些改進能使得該方法更加適用于三維醫(yī)學圖像的分割需求以及實際的臨床應(yīng)用環(huán)境。



圖1為本發(fā)明方法總體流程圖。
圖2為將原始三維醫(yī)學圖像初始化為網(wǎng)絡(luò)流圖的具體流程圖。
圖3為本發(fā)明方法以2x2x2圖像上的分割的實施例示意圖。
圖4為基于CUDA的最大流方法具體實現(xiàn)的流程圖。
圖5為分成兩個階段的relabel操作過程示意圖。

具體實施例方式 本發(fā)明提出的一種基于顯卡并行計算的三維醫(yī)學圖像的快速分割方法，其特征在于，該方法采用基于CUDA技術(shù)加速的Graph-Cuts方法對待分割的原始三維醫(yī)學圖像(從醫(yī)療器械(例如CT機器或是MRI機器)中直接獲得)進行分割處理；如圖1所示，包括以下步驟 1)將待分割的原始三維醫(yī)學圖像初始化為一個網(wǎng)絡(luò)流圖； 2)在網(wǎng)絡(luò)流圖上采用改進的最大流方法進行分割處理，得到三維醫(yī)學圖像的分割結(jié)果； 上述步驟1)將待分割的原始三維醫(yī)學圖像初始化為一個網(wǎng)絡(luò)流圖，如圖2所示，具體包括以下步驟 11)首先根據(jù)原始三維醫(yī)學圖像構(gòu)造出由節(jié)點、端點及有權(quán)邊組成的網(wǎng)絡(luò)流圖的基本架構(gòu)將待分割的原始三維醫(yī)學圖像的所有三維體素都作為網(wǎng)絡(luò)流圖中的節(jié)點，0，1，2，3，…，N-1，N為一個自然數(shù)，表示總共的節(jié)點數(shù)，用邊將相鄰節(jié)點相連接，并為每條邊分配一個容量(即這些有權(quán)邊的權(quán)值，該容量值反映這些節(jié)點對應(yīng)的三維體素值間的梯度信息)，使其成為有權(quán)邊；同時在各節(jié)點的上下兩端添加兩個端點S和T，并將所有節(jié)點分別與兩個端點相連，并為每條邊分配一個容量(即這些有權(quán)邊的權(quán)值，該容量值用于判斷該節(jié)點應(yīng)當在分割結(jié)果中屬于目標器官還是背景)；上述容量值均為非負實數(shù)； 上述步驟11)中網(wǎng)絡(luò)流圖中有權(quán)邊的容量值具體分配如下 定義連接相鄰節(jié)點p和q的有權(quán)邊的容量值為式(1)所示 式(1)中C(p)，C(q)分別是p，q節(jié)點的三維體素值；從式(1)可知，p和q間三維體素之間的差值越大，則該有權(quán)邊的容量值就越小(最終的分割邊界將更傾向于容量值較小的有權(quán)邊)。節(jié)點間有權(quán)邊的流量值與容量值B(p，q)的關(guān)系滿足f(p，q)＝-f(q，p)，|f(p，q)|≤B(p，q)，其中，f(p，q)表示由節(jié)點p流出后流入節(jié)點q的流量，而f(q，p)則表示由節(jié)點q流出后流入節(jié)點p的流量； 定義連接節(jié)點p和端點S和T的有權(quán)邊，其容量值為式(2)所示 R(p，T)＝-lnPr(C(p)|hist″obj″) (2) R(p，S)＝-lnPr(C(p)|hist″bkg″) 上式中hist“obj”和hist“bkg”分別表示預(yù)先獲得的目標器官及背景的三維體素直方圖(例如分別統(tǒng)計手工勾畫或是使用其它分割算法分割出的三維醫(yī)學圖像中目標器官和背景部分中所有三維體素的體素值)，Pr(C(p)|hist“obj”)和Pr(C(p)|hist“bkg”)則分別表示節(jié)點p最終將位于分割目標以及背景中的可能性(體素值與分割目標直方圖接近的節(jié)點就更容易被分割到最終的目標中)；對于節(jié)點與端點間有權(quán)邊的流量值與容量值R(p，T)的關(guān)系滿足|f(S，p)|≤R(p，S)，|f(p，T)|≤R(p，T)，其中，f(S，p)表示從端點S流出后流入節(jié)點p的流量，而f(p，T)表示從節(jié)點p流出后流入端點T的流量。
圖3為以N＝8為例的網(wǎng)絡(luò)流圖基本架構(gòu)，圖中節(jié)點0，1，2，…，7為三維體素對應(yīng)的節(jié)點，上方的S和下方的T為添加的端點。節(jié)點間的有權(quán)邊由無向?qū)嵕€表示，端點和節(jié)點間的有權(quán)邊由有向虛線表示，虛線上箭頭的指向表示在后續(xù)的分割方法中需要找到網(wǎng)絡(luò)流圖中從端點S到端點T的最大流。所有有權(quán)邊的容量值大小由線段的粗細程度來進行示意，盡管圖中的線段粗細程度僅為有限的幾級，但實際計算中容量值根據(jù)上述的等式可能取值為任意大小的非負實數(shù)。
12)根據(jù)網(wǎng)絡(luò)流圖的基本構(gòu)架構(gòu)造網(wǎng)絡(luò)流圖的邏輯結(jié)構(gòu)為每一個節(jié)點和端點分配一個整數(shù)值作為該節(jié)點的標號(不是位置的順序號)，設(shè)初始時端點S的標號值為N，其它所有節(jié)點和端點的標號值為0；同時為每條有權(quán)邊分配一個絕對值不超過該邊容量值的“流量值”以反映該邊的流量，初始時設(shè)所有有權(quán)邊的流量值均為0；即構(gòu)成了網(wǎng)絡(luò)流圖的邏輯結(jié)構(gòu)。
13)構(gòu)造網(wǎng)絡(luò)流圖的存儲結(jié)構(gòu)從而獲得完成實際存儲的網(wǎng)絡(luò)流圖 由于在網(wǎng)絡(luò)流圖中，每個節(jié)點與六個節(jié)點(沿三維空間的三個垂直的坐標軸方向各有兩個相鄰節(jié)點)和兩個端點相鄰，而本發(fā)明應(yīng)用場合中的原始三維醫(yī)學圖像的大小可能非常大(512x512x120)，所以如果直接在每個節(jié)點中存儲與其鄰接的所有邊的索引所需的空間開銷太大。
本發(fā)明構(gòu)造網(wǎng)絡(luò)流圖的存儲結(jié)構(gòu)為基于節(jié)點的由五個三維數(shù)組CoEX，CoEY，CoEZ，CoES和CoET組成，將網(wǎng)絡(luò)流圖中的所有的有權(quán)邊的容量值存儲于五個三維數(shù)組CoEX，CoEY，CoEZ，CoES和CoET中，坐標為(a，b，c)的節(jié)點與鄰接節(jié)點(a-1，b，c)，(a+1，b，c)，(a，b-1，c)，(a，b+1，c)，(a，b，c-1)，(a，b，c+1)以及端點S和T間的有權(quán)邊集即為{CoEX(a，b，c)，CoEX(a+1，b，c)，CoEY(a，b，c)，CoEY(a，b+1，c)，CoEZ(a，b，c)，CoEZ(a，b，c+1)，CoES(a，b，c)，CoET(a，b，c)}，通過節(jié)點在原始三維醫(yī)學圖像中的坐標來獲得與其相連接的有權(quán)邊的索引；網(wǎng)絡(luò)流圖中的流量值也以與上述存儲容量值完全相同的方法存儲于五個三維數(shù)組FoEX，F(xiàn)oEY，F(xiàn)oEZ，F(xiàn)oES和FoET中(即這五個數(shù)組(FoEX等)的存儲方式與上述的五個數(shù)組(CoEX等)相同，即利用坐標來找有權(quán)邊的索引)；這樣就得到了用于實際存儲的網(wǎng)絡(luò)流圖。
上述步驟2)在實際存儲的網(wǎng)絡(luò)流圖上采用改進的最大流方法進行分割處理，如圖4所示，具體包括以下步驟 21)對實際存儲的網(wǎng)絡(luò)流圖初始化活躍節(jié)點作為分割的啟動條件，該活躍節(jié)點為流入該節(jié)點的流量值之和大于流出該節(jié)點的流量值之和的節(jié)點(這部分多出的流量值稱為該節(jié)點的“超額流量”)；令網(wǎng)絡(luò)流圖中所有節(jié)點與端點S間的有權(quán)邊的流量值等于該邊的容量值，即將FoES數(shù)組中的所有元素賦值為CoES中對應(yīng)元素的值，使網(wǎng)絡(luò)流圖中所有的節(jié)點均成為活躍節(jié)點； 22)對實際存儲的網(wǎng)絡(luò)流圖進行基于CUDA實現(xiàn)的并行化改進的(即分為兩個階段并行進行的)relabel操作以更新所有活躍節(jié)點的標號 首先僅對偶節(jié)點并行進行relabel操作(即節(jié)點的三維坐標(a，b，c)滿足a+b+c為偶的節(jié)點稱為“偶節(jié)點”)，如圖5(a)所示，圖中編號為0，3，5，6的節(jié)點三維坐標分別為(0，0，0)，(0，1，1)，(1，0，1)，(1，1，0)，為偶節(jié)點，圖中用有向線段表示當前處理的偶節(jié)點尋找相鄰節(jié)點的方向；而余下的節(jié)點(稱為“奇節(jié)點”)的并行relabel操作將在下一階段進行，如圖5(b)所示，編號1，2，4，7的節(jié)點三維坐標分別為(0，1，0)，(0，0，1)，(1，0，0)，(1，1，1)，為奇節(jié)點，圖中用有向線段表示當前處理的奇節(jié)點尋找相鄰節(jié)點的方向。分別對偶節(jié)點和奇節(jié)點先后進行改進的relabel操作的目的是為了避免對兩個相鄰節(jié)點同時進行relabel操作時可能出現(xiàn)的沖突。舉例而言，考慮只有兩個節(jié)點v和w的圖像，這兩個節(jié)點的初始標號均為0且由一條有權(quán)邊相連。如果同時對這兩個節(jié)點進行relabel操作，則v將被標號為1，因為其標號最小的鄰居為w(標號為0)，而基于同樣的原理w也會被標號為1，這樣在后續(xù)的標號過程中這兩個節(jié)點的標號會一直保持相同的數(shù)值遞增下去，從而導致下一步無法進行任何push操作(這是因為進行push操作的兩個節(jié)點的標號值必須相差1)。使用本發(fā)明給出的機制，沒有任何一對相鄰節(jié)點會被同時進行relabel操作，從而使方法可以順利的運行到后續(xù)步驟。
上述對偶節(jié)點和奇節(jié)點具體進行的改進的relabel操作的步驟完全相同，均為 先找出活躍節(jié)點v的所有鄰接邊的容量值和流量值(假設(shè)v的坐標為(a，b，c)，則其各鄰接邊的容量值分別為CoEX(a，b，c)，CoEX(a+1，b，c)，CoEY(a，b，c)，CoEY(a，b+1，c)，CoEZ(a，b，c)，CoEZ(a，b，c+1)，CoES(a，b，c)，CoET(a，b，c)；各鄰接邊的流量值為FoEX(a，b，c)，F(xiàn)oEX(a+1，b，c)，F(xiàn)oEY(a，b，c)，F(xiàn)oEY(a，b+1，c)，F(xiàn)oEZ(a，b，c)，F(xiàn)oEZ(a，b，c+1)，F(xiàn)oES(a，b，c)，F(xiàn)oET(a，b，c))； 令其中容量值大于流量值的有權(quán)邊所連接的鄰接節(jié)點和端點構(gòu)成集合W，即W＝{w|f(v，w)＜B(v，w)‖f(v，w)＜R(v，w)}，其中f(v，w)表示節(jié)點v和w間有權(quán)邊的流量值，B(v，w)表示w為節(jié)點時v，w間有權(quán)邊的容量值，R(v，w)表示w為端點時v，w間有權(quán)邊的容量值； 令所有活躍節(jié)點按固定順序與W中的點并行進行比較(例如先比較沿x軸的兩個節(jié)點，接下來沿y軸和z軸比較，最后比較端點S和T)；將活躍節(jié)點的標號更新為這些節(jié)點和端點中標號最小者的標號值加一，即d(v)＝min{d(w)+1|w∈W}，其中d(v)表示節(jié)點v的標號； 上述的并行比較具體實現(xiàn)方法為給每個活躍節(jié)點分配一個CUDA線程，然后令多個線程在GPU上同時運行(例如將每512個線程組成一個CUDA架構(gòu)下的block，每512個block組成CUDA架構(gòu)下的一個grid來并行處理)，從而讓多個活躍節(jié)點同時進行relabel操作以減少運行時間； 23)對網(wǎng)絡(luò)流圖中進行基于CUDA實現(xiàn)的兩階段并行化改進的push操作 首先僅對那些偶節(jié)點并行進行push操作，然后對余下的奇節(jié)點并行進行push操作；以將活躍節(jié)點的超額流量轉(zhuǎn)移到一個標號值小一的相鄰節(jié)點或端點上分別對偶節(jié)點和奇節(jié)點先后進行改進的push操作是為了避免同時更新多個節(jié)點超額流量時所可能產(chǎn)生的沖突。舉例而言，如果三個相互鄰接且標號分別為2，1，0的活躍節(jié)點同時進行push操作，標號為2的節(jié)點將會把一定的超額流量e1轉(zhuǎn)移到標號為1的節(jié)點，而此節(jié)點又會將一定的超額流量e2轉(zhuǎn)移到標號為0的節(jié)點。此時中間節(jié)點的超額流量值就可能因為需要同時進行增加e1和減少e2的操作而產(chǎn)生沖突。
上述對偶節(jié)點和奇節(jié)點具體進行的改進的push操作完全相同，具體步驟為 先找出活躍節(jié)點v的所有鄰接邊的容量值和流量值(假設(shè)v的坐標為(a，b，c)，則其各鄰接邊的容量值分別為CoEX(a，b，c)，CoEX(a+1，b，c)，CoEY(a，b，c)，CoEY(a，b+1，c)，CoEZ(a，b，c)，CoEZ(a，b，c+1)，CoES(a，b，c)，CoET(a，b，c)；各鄰接邊的流量值為FoEX(a，b，c)，F(xiàn)oEX(a+1，b，c)，F(xiàn)oEY(a，b，c)，F(xiàn)oEY(a，b+1，c)，F(xiàn)oEZ(a，b，c)，F(xiàn)oEZ(a，b，c+1)，F(xiàn)oES(a，b，c)，F(xiàn)oET(a，b，c))； 令其中與v連接的有權(quán)邊容量值大于流量值，且標號值比v小一的鄰接節(jié)點和端點構(gòu)成集合W，即W＝{w|(f(v，w)＜B(v，w)‖f(v，w)＜R(v，w))∧d(v)＝d(w)+1}，其中f(v，w)表示節(jié)點v和w間有權(quán)邊的流量值，B(v，w)表示w為節(jié)點時v，w間有權(quán)邊的容量值，R(v，w)表示w為端點時v，w間有權(quán)邊的容量值，d(v)表示節(jié)點v的標號； 令所有活躍節(jié)點按固定順序與W中的點并行進行比較(例如先比較沿x軸的兩個節(jié)點，接下來沿y軸和z軸比較，最后比較端點S和T)，然后將超額流量δ＝min(e(v)，B(v，w)-f(v，w)‖R(v，w)-f(v，w))沿該對應(yīng)邊傳遞給集合W中的各節(jié)點和端點，即f(v，w)＝f(v，w)+δ；f(w，v)＝f(w，v)-δ；e(v)＝e(v)-δ；e(w)＝e(w)+δ，其中e(v)表示當前節(jié)點v的超額流量值； 上述的并行比較具體實現(xiàn)方法為給每個活躍節(jié)點分配一個CUDA線程，然后令多個線程同時運行(例如將每512個線程組成一個CUDA架構(gòu)下的block，每512個block組成CUDA架構(gòu)下的一個grid來并行處理)，從而讓多個活躍節(jié)點同時進行push操作以減少運行時間。
24)當活躍節(jié)點的超額流量值被更新為0，則該活躍節(jié)點就變?yōu)槠胀ü?jié)點，檢查網(wǎng)絡(luò)流圖中是否還存在活躍節(jié)點；如果依然存在，則回轉(zhuǎn)到22)步驟繼續(xù)進行relabel和push操作；否則若網(wǎng)絡(luò)流圖中不再存在活躍節(jié)點，網(wǎng)絡(luò)流圖中將得到從S到T的最大流，最大流對應(yīng)的最小割切中容量值等于流量值的邊即構(gòu)成了分割目標器官與背景的邊界，這些邊所連接的三維體素的集合就將原始的三維醫(yī)學圖像分割為了目標器官和背景兩部分，即得到了分割的結(jié)果。
本發(fā)明實現(xiàn)的方法已經(jīng)在不同尺寸下的人類器官的醫(yī)學圖像上進行了方法的測試，表1給出了其中的一些實驗結(jié)果。本發(fā)明使用的CUDA是2.3版本，GPU設(shè)備則是NVIDIA公司的GTX275顯卡，其具有240個流處理單元以及896M字節(jié)的內(nèi)存，而作為比較的Intel雙核CPU主頻是2.50GHz。如表1示 表1.本發(fā)明方法和傳統(tǒng)CPU實現(xiàn)的運行時間開銷的比較 表1所示的結(jié)果表明本發(fā)明可以在10-20倍快于傳統(tǒng)CPU實現(xiàn)的速度下獲得精確的分割結(jié)果。與其它相關(guān)的研究相比本發(fā)明在運行速度方面的優(yōu)勢也是很明顯的。該結(jié)果中對每張原始二維醫(yī)學圖像的平均分割時間遠小于Lee等人的方法中3.35s的開銷。此外，由于網(wǎng)絡(luò)流圖是被存于紋理內(nèi)存中進行處理，本發(fā)明可以利用歸一化的坐標而不是實際的體素坐標來以任意分辨率獲取圖像的信息。也就是說，這里可以構(gòu)造任意大小的網(wǎng)絡(luò)流圖而無需考慮原始三維醫(yī)學圖像的原始尺寸。舉例而言，假設(shè)圖像的原始大小為512x512x36，本發(fā)明可以通過歸一化坐標(0.5，0.4，0.1)來獲取通過插值來得到的原圖位于(512×0.5，512×0.4，36×0.1)處的CT數(shù)值。因此當原始三維醫(yī)學圖像尺寸較大的情況下也能構(gòu)造一個相對小的網(wǎng)絡(luò)流圖。借助這一特性，本發(fā)明可以通過犧牲一些分割的精度來為輸入圖像很大的情形節(jié)省相當?shù)拇鎯π枨蠛瓦\行時間。
權(quán)利要求
1.一種基于顯卡并行計算的三維醫(yī)學圖像的快速分割方法，其特征在于，該方法采用基于CUDA技術(shù)加速的Graph-Cuts方法對待分割的原始三維醫(yī)學圖像進行分割處理；包括以下步驟
1)將待分割的原始三維醫(yī)學圖像初始化為一個網(wǎng)絡(luò)流2)在網(wǎng)絡(luò)流圖上采用改進的最大流方法進行分割處理，得到三維醫(yī)學圖像的分割結(jié)果；
上述步驟1)將待分割的原始三維醫(yī)學圖像初始化為一個網(wǎng)絡(luò)流圖，具體包括以下步驟
11)首先根據(jù)原始三維醫(yī)學圖像構(gòu)造出由節(jié)點、端點及有權(quán)邊組成的網(wǎng)絡(luò)流圖的基本架構(gòu)；
12)根據(jù)網(wǎng)絡(luò)流圖的基本構(gòu)架構(gòu)造網(wǎng)絡(luò)流圖的邏輯結(jié)構(gòu)為每一個節(jié)點和端點分配一個整數(shù)值作為該節(jié)點的標號，設(shè)初始時端點S的標號值為N，其它所有節(jié)點和端點的標號值為0；同時為每條有權(quán)邊分配一個絕對值不超過該邊容量值的流量值，初始時設(shè)所有有權(quán)邊的流量值均為0；即構(gòu)成網(wǎng)絡(luò)流圖的邏輯結(jié)構(gòu)；
13)構(gòu)造網(wǎng)絡(luò)流圖的存儲結(jié)構(gòu)從而獲得完成實際存儲的網(wǎng)絡(luò)流圖
所述存儲結(jié)構(gòu)為基于節(jié)點的由兩組五個三維數(shù)組組成，將網(wǎng)絡(luò)流圖中的所有的有權(quán)邊的容量值存儲于第一組的五個三維數(shù)組CoEX，CoEY，CoEZ，CoES和CoET中，將網(wǎng)絡(luò)流圖中的流量值存儲于第二組的五個三維數(shù)組FoEX，F(xiàn)oEY，F(xiàn)oEZ，F(xiàn)oES和FoET中，通過節(jié)點在原始三維醫(yī)學圖像中的坐標來獲得與其相連接的有權(quán)邊的索引；這樣就得到了用于實際存儲的網(wǎng)絡(luò)流上述步驟2)在實際存儲的網(wǎng)絡(luò)流圖上采用改進的最大流方法進行分割處理，具體包括以下步驟
21)對實際存儲的網(wǎng)絡(luò)流圖初始化活躍節(jié)點作為分割的啟動條件，該活躍節(jié)點為流入該節(jié)點的流量值之和大于流出該節(jié)點的流量值之和的節(jié)點；令網(wǎng)絡(luò)流圖中所有節(jié)點與端點S間的有權(quán)邊的流量值等于該邊的容量值，使網(wǎng)絡(luò)流圖中所有的節(jié)點均成為活躍節(jié)點；
22)對實際存儲的網(wǎng)絡(luò)流圖進行基于CUDA實現(xiàn)的兩階段并行relabel操作
首先僅對偶節(jié)點并行進行relabel操作，偶節(jié)點為節(jié)點三維坐標(a，b，c)滿足a+b+c；然后對余下的奇節(jié)點并行進行relabel操作；以更新所有活躍節(jié)點的標號
23)對完成relabel操作后的網(wǎng)絡(luò)流圖進行基于CUDA實現(xiàn)的兩階段并行push操作首先偶節(jié)點并行進行push操作，然后對余下的奇節(jié)點并行進行push操作；以將活躍節(jié)點的超額流量轉(zhuǎn)移到一個標號值小一的相鄰節(jié)點或端點上
24)當活躍節(jié)點的超額流量值被更新為0，則該活躍節(jié)點就變?yōu)槠胀ü?jié)點，檢查網(wǎng)絡(luò)流圖中是否還存在活躍節(jié)點；如果依然存在，則回轉(zhuǎn)到22)步驟繼續(xù)進行relabel和push操作；否則若網(wǎng)絡(luò)流圖中不再存在活躍節(jié)點，網(wǎng)絡(luò)流圖中將得到從S到T的最大流，最大流對應(yīng)的最小割切中容量值等于流量值的邊即構(gòu)成了分割目標器官與背景的邊界，即得到了分割的結(jié)果。
2.如權(quán)利要求1所述方法，其特征在于，所述步驟11)中網(wǎng)絡(luò)流圖中有權(quán)邊的容量值具體分配如下
定義連接相鄰節(jié)點p和q的有權(quán)邊的容量值為式(1)所示
式(1)中C(p)，C(q)分別是p，q節(jié)點的三維體素值；節(jié)點間有權(quán)邊的流量值與容量值B(P，q)的關(guān)系滿足f(P，q)＝-f(q，P)，|f(p，q)|≤B(P，q)，其中，f(P，q)表示由節(jié)點p流出后流入節(jié)點q的流量，而f(q，P)則表示由節(jié)點q流出后流入節(jié)點p的流量；
定義連接節(jié)點p和端點S和T的有權(quán)邊，其容量值為式(2)所示
R(p，T)＝-InPr(C(p)|hist″obj″)
(2)
R(p，S)＝-InPr(C(p)|hist″bkg″)
上式中hist“obj”和hist“bkg”分別表示預(yù)先獲得的目標器官及背景的三維體素直方圖，Pr(C(p)|hist“obj”)和Pr(C(p)|hist“bkg”)則分別表示節(jié)點p最終將位于分割目標以及背景中的可能性；對于節(jié)點與端點間有權(quán)邊的流量值與容量值R(p，T)的關(guān)系滿足|f(S，p)|≤R(p，S)，|f(p，T)|≤R(p，T)，其中，f(S，p)表示從端點S流出后流入節(jié)點p的流量，而f(p，T)表示從節(jié)點p流出后流入端點T的流量。
3.如權(quán)利要求1所述方法，其特征在于，所述步驟22)中對偶節(jié)點和奇節(jié)點并行進行的relabel操作的步驟完全相同，均為
先找出活躍節(jié)點v的所有鄰接邊的容量值和流量值，設(shè)v的坐標為(a，b，c)，則其各鄰接邊的容量值分別為CoEX(a，b，c)，CoEX(a+1，b，c)，CoEY(a，b，c)，CoEY(a，b+1，c)，CoEZ(a，b，c)，CoEZ(a，b，c+1)，CoES(a，b，c)，CoET(a，b，c)；各鄰接邊的流量值為FoEX(a，b，c)，F(xiàn)oEX(a+1，b，c)，F(xiàn)oEY(a，b，c)，F(xiàn)oEY(a，b+1，c)，F(xiàn)oEZ(a，b，c)，F(xiàn)oEZ(a，b，c+1)，F(xiàn)oES(a，b，c)，F(xiàn)oET(a，b，c))；
令其中容量值大于流量值的有權(quán)邊所連接的鄰接節(jié)點和端點構(gòu)成集合W，即W＝{w|f(v，w)＜B(v，w)||f(v，w)＜R(v，w)}，其中f(v，w)表示節(jié)點v和w間有權(quán)邊的流量值，B(v，w)表示w為節(jié)點時v，w間有權(quán)邊的容量值，R(v，w)表示w為端點時v，w間有權(quán)邊的容量值；
令所有活躍節(jié)點按固定順序與W中的點并行進行比較；將活躍節(jié)點的標號更新為這些節(jié)點和端點中標號最小者的標號值加一，即d(v)＝min{d(w)+1|w∈W}，其中d(v)表示節(jié)點v的標號；
所述的并行比較具體實現(xiàn)方法為給每個活躍節(jié)點分配一個CUDA線程，然后令多個線程在GPU上同時運行，使多個活躍節(jié)點同時進行relabel操作，以減少運行時間；
4.如權(quán)利要求1所述方法，其特征在于，所述步驟23)對偶節(jié)點和奇節(jié)點并行進行push操作完全相同，均為
先找出活躍節(jié)點v的所有鄰接邊的容量值和流量值，設(shè)v的坐標為(a，b，c)，則其各鄰接邊的容量值分別為CoEX(a，b，c)，CoEX(a+1，b，c)，CoEY(a，b，c)，CoEY(a，b+1，c)，CoEZ(a，b，c)，CoEZ(a，b，c+1)，CoES(a，b，c)，CoET(a，b，c)；各鄰接邊的流量值為FoEX(a，b，c)，F(xiàn)oEX(a+1，b，c)，F(xiàn)oEY(a，b，c)，F(xiàn)oEY(a，b+1，c)，F(xiàn)oEZ(a，b，c)，F(xiàn)oEZ(a，b，c+1)，F(xiàn)oES(a，b，c)，F(xiàn)oET(a，b，c))；
令其中與v連接的有權(quán)邊容量值大于流量值，且標號值比v小一的鄰接節(jié)點和端點構(gòu)成集合W，即W＝{w|(f(v，w)＜B(v，w)||f(v，w)＜R(v，w))∧d(v)＝d(w)+1}，其中f(v，w)表示節(jié)點v和w間有權(quán)邊的流量值，B(v，w)表示w為節(jié)點時v，w間有權(quán)邊的容量值，R(v，w)表示w為端點時v，w間有權(quán)邊的容量值，d(v)表示節(jié)點v的標號；
令所有活躍節(jié)點按固定順序與W中的點并行進行比較，然后將超額流量δ＝min(e(v)，B(v，w)-f(v，w)||R(v，w)-f(v，w))沿該對應(yīng)邊傳遞給集合W中的各節(jié)點和端點，即f(v，w)＝f(v，w)+δ；f(w，v)＝f(w，v)-δ；e(v)＝e(v)-δ；e(w)＝e(w)+δ，其中e(v)表示當前節(jié)點v的超額流量值；
所述的并行比較具體實現(xiàn)方法為給每個活躍節(jié)點分配一個CUDA線程，然后令多個線程同時運行，使多個活躍節(jié)點同時進行push操作，以減少運行時間。
全文摘要
本發(fā)明涉及一種基于顯卡并行計算的三維醫(yī)學圖像的快速分割方法，屬于醫(yī)學圖像處分割領(lǐng)域，該方法首先根據(jù)原始三維醫(yī)學圖像構(gòu)造出由節(jié)點、端點及有權(quán)邊組成的網(wǎng)絡(luò)流圖的基本架構(gòu)；再構(gòu)造網(wǎng)絡(luò)流圖的邏輯結(jié)構(gòu)構(gòu)造網(wǎng)絡(luò)流圖的存儲結(jié)構(gòu)從而獲得完成實際存儲的網(wǎng)絡(luò)流圖對實際存儲的網(wǎng)絡(luò)流圖進行基于CUDA實現(xiàn)的兩階段并行relabel操作和push操作當活躍節(jié)點的超額流量值被更新為0，則該活躍節(jié)點就變?yōu)槠胀ü?jié)點，則網(wǎng)絡(luò)流圖中將得到從S到T的最大流，最大流對應(yīng)的最小割切中容量值等于流量值的邊即構(gòu)成了分割目標器官與背景的邊界，即得到了分割的結(jié)果。本發(fā)明方法能夠在保持高處理精度的前提下大幅提高分割方法的運行速度。
文檔編號G06T7/00GK101783028SQ201010115248
公開日2010年7月21日 申請日期2010年2月26日 優(yōu)先權(quán)日2010年2月26日
發(fā)明者王宏, 楊凡, 翟偉明, 宋亦旭, 趙雁南, 賈培發(fā) 申請人:清華大學