專利名稱：基于gpu加速的三維醫(yī)學圖像顯示方法
技術領域：
本發(fā)明屬于醫(yī)學圖像處理技術領域，土要利用GPU (顯示處理單元)的強大的并行流處 理能力，來對體繪制進行加速，相較亍傳統(tǒng)CPU的體繪制技術，該方法使得體繪制能交互式 實時顯示。
背景技術：
醫(yī)學圖像在醫(yī)生診斷中的輔助作用越來越明顯，但是僅可以從二維截面方向對人體進行 觀察。為了提高醫(yī)療診斷和治療規(guī)劃的準確性和科學性，需要由二維斷層圖像序列轉變?yōu)榫?有直觀立體效果的三維圖像。然而目前的三維醫(yī)療輔助t會斷系統(tǒng)都需要工作站級別的運行平
臺，才能基本滿足實時性的要求，由于價格問題使得這種系統(tǒng)難以推廣。目能國際上計算機
圖形圖像學的科學家們提出將大規(guī)模的數據處理、運算等任務放在GPU運行作為甜沿的研 究方向，具有CPU不可比擬的速度優(yōu)勢，在消費級的GPU硬件上實現大規(guī)模數據的交互式 可視化。
當甜通用的體繪制技術， 一般是基于傳統(tǒng)的CPU進行計算，對于傳統(tǒng)消費級別大眾PC，
很難實時進行交互式顯示。如果要交互式實時顯示運算， 一般要在并行圖形工作站上進行運 行，成本相較大大提高。
發(fā)明內容
本發(fā)明針對現有基于CPU的醫(yī)學圖像顯示方法存在運算速度慢、無法交互式執(zhí)行的弊 端，提供一種基于GPU加速的三維醫(yī)學圖像顯示方法，利用GPU強大的流運算能力進行加 速，從而實現交互式體繪制的效果。
本發(fā)明技術方案如下
基于GPU加速的三維醫(yī)學圖像顯示方法，如圖1所示，包括以下,:p驟 歩驟1:讀入醫(yī)學DICOM圖像序列文件并以體數據的方式保存到系統(tǒng)內存。 體數據是由有順序的2維醫(yī)學DICOM圖像序列構成。將這些醫(yī)學圖像的圖像相關信息 讀入系統(tǒng)內存(圖像分辨率，層間距，圖像像素信息)。為了更好的視覺效果，體數據應該首 先進行預處理，比如用圖像濾波器進行除噪，并且通過插值層間距數據來得到更加細致的效 果。
歩驟2:利用OpenGL或者DirectX的三維圖形庫編程擴展接口函數API,將歩驟1保存在系統(tǒng)內存中的體數據加載入GPU顯存，成為GPU可以訪問的3D紋理。具體接口函數API 可以采用OpenGL庫的擴展函數glTexImage3DEXT來將體數據加載為3D紋理。
歩驟3:計算生成代理幾何體，具體分為以下幾個歩驟
歩驟3-1:計算體數據包絡盒的頂點在視點坐標系下的坐標。
這個操作涉及到一個頂點坐標變換操作，其實質就是將體數據包絡盒的頂點坐標由體數 據的局部坐標系轉換到視點坐標系。首先將體數據包絡盒八個頂點的局部坐標轉換成世界坐 標，然后通過視點世界坐標參數矩陣，將體數據包絡盒八個頂點的世界坐標變換到視點坐標 系中。
歩驟3-2:計算出代理兒何體所包含的多邊形切片數目。
在視點坐標系下，首先計算體數掂包絡盒的八個頂點在視線方向的最大值zma,和最小值
—7m。；然后確定采樣率A力，即相鄰兩層多邊形切片W的距離；最后計算代理幾何體所包含的 多邊形切片數目A^"、為
步驟3-3:計算代理幾何體每張多邊形切片的頂點的世界坐標。
對于代理兒何體中的每張多邊形切片，通過計算每張多邊形切片所在平面與體數據包絡
盒在世界坐標系下的交點，然后將這些交點的出:界坐標和3D紋理坐標順序地(比如沿順時
針方向或逆時針方向)組成每張多邊形切片的頂點數組。其中，計算交點的具體方法是
設體數據在局部坐標系下數據區(qū)間為[-x,,,m, xm, ], [imin, }mll，]，卜z,，, 2mm ]，視點位置為
'視線方向向量為^ : fe,，>W，^),采樣率為厶A 。對于每張多邊形切
片，即垂直于視線方向、間隔為A/7的平行平面與立方體邊界線段切割的交集區(qū)域。 一張多邊
形切片所在平面可以表示成方程
d/廠p0) = 0
根據視線方向和所經過的體數據內點A,的位置，順序計算體數據邊界線段與平面
= 0的交點。由于我們將體數據局部坐標系和世界坐標系重合，極大方便我們的計 算，因為邊界線段都是與坐標軸平行的。
^驟4:對歩驟3所得的每張多邊形切片，逐像素進行光照計算和顏色計算，具體分為
以下幾個歩驟步驟4-1:計算代理幾何中每張多邊形切片的每一像素點的體數據梯度。
—個三維數據/Oc,少,Z)的梯度由一個偏微分梯度算子表示為V/(X,y,Z)—i，！,！)；
&
具體采用基于6鄰域采樣差分的梯度計算方法，計算位于(X,)'，Z)處的梯度V/(X,，^，^):
V/(乂,,7,,zJ =1 /2(,(x,+1'少,，z》—,(x, i,-"'4)，、V''z*)-/化')',(x''X'zw)—/化'少"z",》
歩驟4-2:采用顏色傳遞凼數計算代理兒何中每張多邊形切片的每一像素點的顏色。
由于體數據圖像為灰度圖像，在視覺上效果以及色度分辨率不利于對細節(jié)進行分辨，所
以在本方法中通過顏色傳遞函數進行體數據體素灰度到RGB顏色空間的映射。
顏色傳遞函數是一個類似c二/0)的一維函數，c為輸出的顏色值，v為體數據的灰度值。
顏色傳遞函數的作用是用來將灰度值的體數據映射到RGB顏色空間的體數據?？梢酝ㄟ^一個顏色映射表來進行實現，在GPU里面，可以通過l維紋理來存放這個顏色映射表。在醫(yī)學圖像序列畢.面，灰度數據是從O到N離散存放的。在這早，N是一個山顏色存儲位數i來決定
的自然數# = 2'。通過預計算出函數映射表的1維紋理，通過對這個紋理進行采樣，這樣很容易得到輸出紋理。這段代碼在GPU eg語舎単面只是一個簡單的指令texld(float x)。
步驟4-3:根據Phong光照模型計算代理幾何體中每張多邊形切片的每一像素點的光柵化像素著色值。
反射光分為三個部分，環(huán)境光，漫反射光，鏡面反射光。
環(huán)境光本質是指光源經過環(huán)境等作用，然后間接對物體的影響，其本質是在物體和環(huán)境之間多次反射，最終達到平衡時的一種光強。本質上說，這種光強在物體上各點是不相同的，但是在簡巾.光照模型中近似地認為同一環(huán)境下的環(huán)境光，其光強分布是均勻的，而且環(huán)境光在任何一個方向上的光強強度都相同。表現出來是一個常數。
漫反射光漫反射光的定義如下，假設當光源來自一個方向時，漫反射光將光均勻向各
方向反射并傳播丌米，與反射光的強度與視點無關，因為漫反射光本質是由反射表面的微表
面凹凸不平所引起的，因而漫反射光的空間分布是均勻的。運算的式子如下；假設入射光的強度為/,,，物體表面微元上點P的法向為W，從點P指向入射光的光源的單位向量為Z，
兩者間的夾角為6 ，由光強能量的余弦定律可以計算得到漫反射光強為乙=/,，>< & xcos(S), 0e(O,"/2)在上式子中，《,是與物體有關的光的漫反射系數，系數滿
足0<《，<1 。當丄、V為單位向量時，也可用如下向量運算的形式進行表達
6乙=;x & x &' AO在有多個光源的情況下'可以有如下的表示/rf = &Z ; x (丄,.W)
鏡面反射光最典型的例子就是對于理想鏡面，鏡面的反射光集中在一個方向，并遵守反射定律。將其推廣，對一般的光滑表面，反射光集中在一個空間區(qū)間范圍內，且由反射定律決定的反射方向光強最大。因此，對于同一點來說，從不同位置所觀察到的鏡面反射光強是不同的，這個與視線與鏡面反射光方向有關。鏡面反射光強的式子可表示為
=x《、x cos"W), P e (O,;r / 2)其中X、是與物體表面光學性質有關的鏡面反射系數，P為
視線方向F與反射方向A的夾角，"為反射指數，其本質實質上是反映了物體表面的光澤程度， 一般為1 2000，通過式子可以看出，數目越大物體表面越光滑。鏡面反射光將會在反射方向附近形成很亮的光斑，稱為高光現象。同樣，如果將F和i 都格式化為單位向量，那
么鏡面反射光強可表示為/、 /p.《(i .r)"其中，A 口J由向量反射計算式子i^2W(7V.丄)-丄計算。與甜面漫反射光相向的，當對十多個光源的情形，鏡面反射光強可以表示為
/、=《^[/,(7 ,J/)"]
歩驟5:通過Alpha混合將代理幾何體中所有多邊形切片合成二維醫(yī)學圖像。
在關閉深度測試的情況下，打丌Alpha混合，將將代理幾何體中所有多邊形切片合成三維醫(yī)學圖像。為了提高渲染速度，可以將代理幾何體保存到顯示列表，這樣就能用單--指令渲染多個頂點，從而提高繪制效率。
Alpha混合是一個將代理幾何體若十多邊形層的顏色進行混合的方法。圖形GPU硬件用的光照模型是針對多邊形頂點的，而不能被用在體繪制上面，所以必須關閉光照。并且，由于打丌了 Alpha混合，就必須渲染代理幾何體単面的每一個多邊形層，所以必須關閉深度測試。下面進行簡單解釋。
Alpha混合OpenCH,中的絕大多數特效都與某些類型的(一般是色彩)混合有關?；焐亩x為，將某個象素的顏色和已繪制在屏幕上與其對應的象素顏色相互結合。至于如何結合這兩個顏色則依賴于顏色的alpha通道的分量值，以及/或者所使用的混色函數。Alpha通常是位于顏色值末尾的第4個顏色組成分量。大多數情況都認為Alpha分量代表材料的透明度。這就是說，alpha值為0.0時所代表的材料是完全透明的。alpha值為1.0時所代表的材料則是完全不透明的。其實混合的基本原理是就將要分色的圖像各象素的顏色以及背景顏色均按照RGB規(guī)則各自分離之后，根據 一 圖像的RGB顏色分量+alpha值+背景的RGB顏色分量^l-alpha值) 一 這樣一個簡單公式來混合之后，最后將混合得到的RGB分量重新合并，公式如下Co/or/ma, = Co/or/ra , xa + Co/o;ti x (1 -a) 。 OpenGL按照卜.面的公式計算這兩個象素
的混色結果。這個公式一般會生成混合的透明/半透明的效果。
深度測試像素片斷測試其實就是測試每一個像素，只有通過測試的像素才會被繪制，沒有通過測試的像素則不進行繪制。OpenGL提供了多種測試操作，利用這些操作可以實現一些特殊的效果。"深度測試"的概念在繪制三維場景的時候特別有用。在不使用深度測試的時候，如果我們先繪制一個距離較近的物體，再繪制距離較遠的物體，則距離遠的物體因為后繪制，會把距離近的物體覆蓋掉，這樣的效果并不是我們所希望的。如果使用了深度測試，則情況就會有所不同每當一個像素被繪制，OpenGL就記錄這個像素的"深度"(深度可以理解為該像素距離觀察者的距離。深度值越大，表示距離越遠)，如果有新的像素即將覆蓋原來的像素時，深度測試會檢查新的深度是否會比原來的深度值小。如果是，則覆蓋像素，繪制成功；如果不是，則不會覆蓋原來的像素，繪制被取消。這樣一來，即使我們先繪制比較近的物體，再繪制比較遠的物體，則遠的物體也不會覆蓋近的物體了。實際上，只要存在深度緩沖區(qū)，無論是否啟用深度測試，OpenGL在像素被繪制時都會嘗試將深度數據寫入到緩沖區(qū)內，除非調用了 glDepthMask(GL—FALSE)來禁止寫入。這些深度數據除了用于常規(guī)的測試外，還可以有一些有趣的用途，比如繪制陰影等等。除了深度測試，OpenGL還提供了剪戰(zhàn)測試、Alpha測試和模板測試。
本發(fā)明提供了--種基于GPU加速的醫(yī)學圖像顯示方法，相比于現有的基于CPU的醫(yī)學圖像顯示方法，本發(fā)明具有很高的運算速度，可在普通消費級別大眾PC上實現實時交互式顯示；而無須使用圖形工作站，使得成本大大降低。
圖1為木發(fā)明流程示意圖。
具體實施例方式
本發(fā)明采用上述技術方案，分別針對MRI和CT醫(yī)學DICOM圖像序列進行顯示驗證，均具有較好的顯示效果。需要說明的是，由于MRI圖像山于噪聲較大，需要進行除噪預處理;而CT圖像較為清晰，無需進行除噪預處理，可直接以體數據的方式進行顯示。
權利要求
1、基于GPU加速的三維醫(yī)學圖像顯示方法，包括以下步驟步驟1讀入醫(yī)學DICOM圖像序列文件并以體數據的方式保存到系統(tǒng)內存；步驟2利用OpenGL或者DirectX的三維圖形庫編程擴展接口函數API，將步驟1保存在系統(tǒng)內存中的體數據加載入GPU顯存，成為GPU可以訪問的3D紋理；步驟3計算生成代理幾何體，具體分為以下幾個步驟步驟3-1計算體數據包絡盒的頂點在視點坐標系下的坐標；首先將體數據包絡盒八個頂點的局部坐標轉換成世界坐標，然后通過視點世界坐標參數矩陣，將體數據包絡盒八個頂點的世界坐標變換到視點坐標系中；步驟3-2計算代理幾何體所包含的多邊形切片數目；在視點坐標系下，首先計算體數據包絡盒的八個頂點在視線方向的最大值zmax和最小值zmin；然后確定采樣率Δh，即相鄰兩層多邊形切片間的距離；最后計算代理幾何體所包含的多邊形切片數目NShces為<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><msub> <mi>N</mi> <mi>Shces</mi></msub><mo>=</mo><mfrac> <mrow><mo>(</mo><msub> <mi>z</mi> <mi>max</mi></msub><mo>-</mo><msub> <mi>z</mi> <mi>min</mi></msub><mo>)</mo> </mrow> <mi>Δh</mi></mfrac> </mrow>]]></math></maths>步驟3-3計算代理幾何體每張多邊形切片的頂點的世界坐標；對于代理幾何體中的每張多邊形切片，通過計算每張多邊形切片所在平面與體數據包絡盒在世界坐標系下的交點，然后將這些交點的世界坐標和3D紋理坐標順序地組成每張多邊形切片的頂點數組；步驟4對步驟3所得的每張多邊形切片，逐像素進行光照計算和顏色計算，具體分為以下幾個步驟步驟4-1計算代理幾何中每張多邊形切片的每一像素點的體數據梯度；一個三維數據f(x，y，z)的梯度由一個偏微分梯度算子表示為<maths id="math0002" num="0002" ><math><![CDATA[ <mrow><mo>▿</mo><mi>f</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo></mrow><mo>=</mo><mrow> <mo>(</mo> <mfrac><mrow> <mo>∂</mo> <mi>f</mi></mrow><mrow> <mo>∂</mo> <mi>x</mi></mrow> </mfrac> <mo>,</mo> <mfrac><mrow> <mo>∂</mo> <mi>f</mi></mrow><mrow> <mo>∂</mo> <mi>y</mi></mrow> </mfrac> <mo>,</mo> <mfrac><mrow> <mo>∂</mo> <mi>f</mi></mrow><mrow> <mo>∂</mo> <mi>z</mi></mrow> </mfrac> <mo>)</mo></mrow><mo>;</mo> </mrow>]]></math> id="icf0002" file="A2009100598640002C2.tif" wi="45" he="9" top= "204" left = "139" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>具體采用基于6鄰域采樣差分的梯度計算方法，計算位于(x，y，z)處的梯度 id="icf0003" file="A2009100598640002C3.tif" wi="23" he="4" top= "220" left = "157" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/><maths id="math0003" num="0003" ><math><![CDATA[ <mrow><mo>▿</mo><mi>f</mi><mrow> <mo>(</mo> <msub><mi>x</mi><mi>i</mi> </msub> <mo>,</mo> <msub><mi>y</mi><mi>j</mi> </msub> <mo>,</mo> <msub><mi>z</mi><mi>k</mi> </msub> <mo>)</mo></mrow><mo>=</mo><mn>1</mn><mo>/</mo><mn>2</mn><mrow> <mo>(</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mrow><mi>i</mi><mo>+</mo><mn>1</mn> </mrow></msub><mo>,</mo><msub> <mi>y</mi> <mi>j</mi></msub><mo>,</mo><msub> <mi>z</mi> <mi>k</mi></msub><mo>)</mo> </mrow> <mo>-</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mrow><mi>i</mi><mo>-</mo><mn>1</mn> </mrow></msub><mo>,</mo><msub> <mi>y</mi> <mi>j</mi></msub><mo>,</mo><msub> <mi>z</mi> <mi>k</mi></msub><mo>)</mo> </mrow> <mo>,</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>,</mo><msub> <mi>y</mi> <mrow><mi>j</mi><mo>+</mo><mn>1</mn> </mrow></msub><mo>,</mo><msub> <mi>z</mi> <mi>k</mi></msub><mo>)</mo> </mrow> <mo>-</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>,</mo><msub> <mi>y</mi> <mrow><mi>j</mi><mo>-</mo><mn>1</mn> </mrow></msub><mo>,</mo><msub> <mi>z</mi> <mi>k</mi></msub><mo>)</mo> </mrow> <mo>,</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>,</mo><msub> <mi>y</mi> <mi>j</mi></msub><mo>,</mo><msub> <mi>z</mi> <mrow><mi>k</mi><mo>+</mo><mn>1</mn> </mrow></msub><mo>)</mo> </mrow> <mo>-</mo> <mi>f</mi> <mrow><mo>(</mo><msub> <mi>x</mi> <mi>i</mi></msub><mo>,</mo><msub> <mi>y</mi> <mi>j</mi></msub><mo>,</mo><msub> <mi>z</mi> <mrow><mi>k</mi><mo>-</mo><mn>1</mn> </mrow></msub><mo>)</mo> </mrow> <mo>)</mo></mrow> </mrow>]]></math></maths>步驟4-2采用顏色傳遞函數計算代理幾何中每張多邊形切片的每一像素點的顏色；其中，所述顏色傳遞函數是一個類似c＝f(v)的一維函數，c為輸出的顏色值，v為體數據的灰度值；顏色傳遞函數的作用是用來將灰度值的體數據映射到RGB顏色空間的體數據；步驟4-3根據Phong光照模型計算代理幾何體中每張多邊形切片的每一像素點的光柵化像素著色值；所述Phong光照模型中反射光分為三個部分，環(huán)境光，漫反射光，鏡面反射光；所述環(huán)境光采用同一環(huán)境下的環(huán)境光，其光強分布是均勻的，而且環(huán)境光在任何一個方向上的光強強度都相同；所述漫反射光的光強為Id＝Ip×Kd×cos(θ)，θ∈(0，π/2)；其中Ip為入射光的強度；Kd是與物體有關的光的漫反射系數，且滿足0＜Kd＜1；θ為物體表面微元上點P的法向N與點P指向入射光的光源的單位向量L的夾角；所述鏡面反射光的光強為Id＝Ip×Ks×cosn(θ′)，θ′∈(0，π/2)；其中Ip為入射光的強度；Ks是與物體表面光學性質有關的鏡面反射系數，θ′為視線方向V與反射方向R的夾角，n為反射指數；步驟5通過Alpha混合將代理幾何體中所有多邊形切片合成三維醫(yī)學圖像；在關閉深度測試的情況下，打開Alpha混合，將將代理幾何體中所有多邊形切片合成三維醫(yī)學圖像。為了提高渲染速度，可以將代理幾何體保存到顯示列表，這樣就能用單一指令渲染多個頂點，從而提高繪制效率。
2、 根據權利要求1所述的基于GPU加速的二維醫(yī)學圖像顯示方法，其特征在于，歩驟 1讀入醫(yī)學D1C0M圖像序列文件并以體數據的方式保存到系統(tǒng)內存時，為了更好的視覺效 果，先對體數據應改進行除U呆或插值層間數據的預處理。
3、 根據權利要求1所述的基于GPU加速的三維醫(yī)學圖像顯示方法，其特征在于，歩驟 2中具體接口函數API采用OpenGL庫的擴展函數glTexImage3DEXT來將體數據加載為3D 紋理。
4、 根據權利要求1所述的基于GPU加速的三維醫(yī)學圖像顯示方法，其特征在于，歩驟 3-3中交點的具體計算方法是設體數據在局部坐標系下數據區(qū)間為[-xmil', x■ ]， [-ym,n. ym,n ]. [-—7,細,zmm ],視點位置為^^:",,<,,，_y_，^. .)，視線方向向量為^:(.、,，U,〃)，采樣率為A/7; —張多邊形切片所 在平面可以表示成方程 7r .(x — p0) = 0根據視線方向和所經過的體數據內點/;。的位置，順序計算體數據邊界線段與平面 (A'r (x — / 0) = 0白勺交點。
全文摘要
基于GPU加速的三維醫(yī)學圖像顯示方法，屬于醫(yī)學圖像處理技術領域。首先將醫(yī)學DICOM圖像序列文件以體數據的方式保存到系統(tǒng)內存；然后利用OpenGL或者DirectX的三維圖形庫編程擴展接口函數API，將體數據加載入GPU顯存；再計算生成代理幾何體，并代理幾何體中的多邊形切片逐像素進行光照計算和顏色計算；最后通過Alpha混合將代理幾何體中所有多邊形切片合成三維醫(yī)學圖像。本發(fā)明相比于現有的基于CPU的醫(yī)學圖像顯示方法，本發(fā)明具有很高的運算速度，可在普通消費級別大眾PC上實現實時交互式顯示；而無須使用圖形工作站，使得成本大大降低。
文檔編號G06T15/00GK101593345SQ20091005986
公開日2009年12月2日 申請日期2009年7月1日 優(yōu)先權日2009年7月1日
發(fā)明者帆 張, 梅 解 申請人:電子科技大學