一種基于渲染的快速體數據骨架提取方法
【專利摘要】本發明涉及一種基于渲染的快速體數據骨架提取方法，包括以下步驟：在上下前后左右6個視點下，用RayCasting直接體繪制算法對三維體數據進行渲染，得到6個2D圖像；對所有直接體繪制獲得的2D圖像提取其輪廓及部分關鍵特征線；將獲得的2D輪廓線及關鍵特征線反投影回3D空間，獲得3D候選曲線；對空間采樣，計算空間采樣點的曲線密度，構建曲線密度場；獲得曲線密度場后，用三維Sobel算子采樣點的梯度；根據采樣點的梯度值，用最小二乘方法計算采樣點的曲線方向；將采樣點的曲線方向賦值給其所對應的三維體數據的體素；根據體素的方向，構建包圍盒，計算包圍盒體素平均值獲得表征節點；連接表征節點并平滑曲線最終得到三維體數據骨架。
【專利說明】
【技術領域】
[0001] 本發明屬于骨架提取及體繪制【技術領域】，具體涉及基于渲染的快速體數據骨架提 取方法。 一種基于渲染的快速體數據骨架提取方法

【背景技術】
[0002] 骨架作為形狀表示的一種有效形式在三維模型的各個研究領域中被廣泛采用。 Bluml967年給出了骨架的最初定義：骨架（中軸）是模型內部各個最大內切球中心的集 合。它還有一個grassfire的模擬定義，從模型表面開始點火，各個方向上的火的相遇點所 構成的集合。因為模型的骨架很好的保留了模型的拓撲連接性及其幾何形狀信息，所以經 常被用于碰撞檢測、三維動畫、模型渲染、模型表面重建、模型檢索、三維形狀識別與檢索、 模型編輯、簡化以及三維分割等應用中，也有研究人員采用骨架為模型的分解做矯正。不同 的應用，對于骨架應該保存的信息要求不完全相同，故而抽取思路也不完全相同。
[0003] 對三維模型骨架的研究由來已久，出現過很多方法，有的是源于對二維圖像的擴 展，有的是針對三維模型提出的，大體上來說，面向三維網格模型或者三維體素模型的骨架 提取方法主要包括：細化、距離變換、分解等。
[0004] 提取三維體素模型骨架的拓撲細化方法的基本思想是：由外向內，層層剝離不改 變模型拓撲特性的體素來實現骨架抽取，它檢測物體數據的拓撲相關性，而不是形狀的幾 何特征。這種方法在保持物體拓撲性質的同時，在離散域中反復的逐層剝離去物體的表面， 其中判斷一個體素是否需要剝離的檢測過程比較耗時，這種發方法需要完全分割好的、緊 湊的并且相互連接的物體模型。Bertrand提出一種簡單的并行細化方法抽取骨架。Ju采 用邊細化邊去除冗余體素的方法抽取骨架。Melada提出結合歐式距離測度與細化過程的骨 架抽取方法。但這些方法產生的骨架包含雜枝，甚至破壞了骨架的連續性。
[0005] 基于距離變換的骨架提取方法：它一般的計算對象要求也必須是體素表示的模 型，通過計算每個體素的距離來求取模型的脊點。Zhou提出了邊界距離變換和種子點距離 變換，采用體素編碼的方式抽取并優化骨架結構，但骨架在分支處的中心性較差。Wan的方 法采用精確的邊界距離變換，建立極小擴張樹提取骨架，從而使得到的骨架點盡量不偏離 中心，但會出現中心線的"擾動"現象。基于Level Set的骨架提取方法穩定性較好，也具有 較高的拓撲無關性，能夠有效消除尖點、骨架斷裂等問題，但是即便采用了 Fast Marching Method策略，其計算復雜度依然很高。Fan引入snake模型調整骨架點的位置，改善了骨架 點的中心性，但復雜度也同時增大。
[0006] 相對于基于細化和距離變換兩種方法而言，基于形狀分解或者分割的三維骨架提 取方法將骨架提取和形狀分割結合起來，骨架提取與形狀分解二者相輔相成具有顯著優 勢。基于拓撲連接信息，先行對三維模型進行分解；基于分解結果可以降低骨架提取的復雜 度、提高骨架的精確度。按照拓撲結構，正確分解三維模型能夠保證快速準確地獲取骨架， 反之，骨架作為拓撲結構的參考能夠保證模型分解的準確性。基于分解的方法源于對形體 的認識，Hoffman和Richards把形體自然分解成基本部件，把分解用于分類。Rom等提出利 用中軸逐步分解二維形體的方法，先分解形體，在提取分解區域的中軸，而中軸是二維形體 的自然表示，但是這些方法難以提取中心良好的三維立體形體的分級骨架。Minimal Rule Lien等提出基于模型近似圖面體分解，通過計算三維模型頂點的凹凸性展開分解，并通過 迭代生成骨架，但得到的骨架中心性較差，也難以建立層次關系。
[0007] 基于Reeb圖思想。該類算法首先在模型上定義一個連續函數，計算每個模型頂點 的函數值，將具有相同函數值的頂點聚合成一個頂點，得到模型的骨架。Hilaga等人提出了 基于Reeb圖思想的多分辨率Reeb圖（MRG)，將Reeb圖用于骨架提取，并將映射函數定義為 頂點到整個模型表面所有頂點的最短距離與面積之乘積的和。Julien Tierny等人使用了 更加合理的標量函數計算方法。
[0008] 亦有學者將二維圖像領域中的Voronoi圖技術引入到三維骨架抽取中，Dey等人 提出過一種利用Voronoi圖直接近似中軸的算法。因為對三維模型而言只有部分Voronoi 頂點能夠匯聚成骨架，所以他通過定義角度和比例這兩個與大小、比重都無關的篩選標準 來實現它的中軸近似算法，并證明該方法能夠保證收斂。Voronoi圖方法對于復雜的三維模 型，理論上可以實現，但計算量龐大。Voronoi圖方法多用于生成多尺度骨架，但對于邊界噪 聲非常敏感，導致Voronoi圖過于稠密，需要額外進行剪枝處理，而且規則復雜，導致該方 法實用性不光。
[0009] 體繪制技術是一種重要的科學可視化技術，采用該技術，不經過復雜的特征算法， 就能夠交互地瀏覽數據內部結構特征。該項技術廣泛應用于醫學、氣象和地質等領域。目前 直接體繪制算法主要有Ray Casting算法、Shear-Warp算法、足跡算法、Texture Mapping 算法等。LevoyM等將體數據遞歸均勻子分成一個八叉樹結構，在光線投射算法中遍歷八叉 樹來跳過包含當前采樣點的最大空區域，以此達到加速的效果。
[0010] 綜上所述，現有的技術存在著計算復雜度大，計算時間長，骨架連續性、中心性差、 有擾動，骨架對噪聲敏感等缺點，本發明提供一種快速、精確、噪聲擾動不敏感的骨架提取 方法。


【發明內容】

[0011] 本發明要解決的技術問題是：克服了現有技術計算時間長、骨架中心性差、對噪聲 敏感的不足，提供一種提取體數據骨架的方法，能夠快速地提取出三維體數據模型的骨架， 并且能得到較好的效果。
[0012] 本發明解決上述技術問題的技術方案為：基于渲染的快速體數據骨架提取方法， 利用NVIDIA GPU硬件的計算能力，提高了計算效率，該方法主要利用具有視覺意義的2D輪 廓線及部分關鍵特征線，并根據能量最小化函數，計算對應的三維特征線，對空間采樣后構 建曲線密度場進而求取采樣點的曲線方向，并將采樣點方向賦值給相應體素以構建包圍盒 得到表征節點，連接這些表征節點并平滑曲線得到體數據骨架。
[0013] 如圖1所示，具體實現如下步驟：
[0014] 1) 3D候選曲線的提取
[0015] 步驟（1)、在X、Y、Z3個坐標軸的正負方向（上下前后左右）上設置六個視點，用 RayCasting直接體繪制算法對三維體數據進行渲染，得到6幅2D圖像；該方法是對圖像的 每一個像素，沿固定方向（通常是視線方向）發射一條光線，光線穿越體數據的整個圖像序 列，并在這個過程中，對圖像序列進行采樣獲取顏色信息，同時依據光線吸收模型將顏色值 進行累加，直至光線穿越整個圖像序列，最后得到的顏色值就是渲染圖像的顏色。該過程可 以在GPU上實現，利用CUDA(計算統一設備架構）的硬件加速，大大提高渲染速度。
[0016] 步驟（2)、對步驟⑴得到的6幅2D圖像，使用輪廓檢測算法提取2D圖像的輪廓， 直接提取的輪廓需要收縮以獲得有效輪廓線，另外為表示更多細節需要使用Canny算子對 2D圖像進行部分關鍵特征線的提取。Canny算子是一階算子，其邊緣檢測可以分為三個步 驟：①高斯平滑函數。目的是為了平滑以消除噪聲；②一階差分卷積模板。目的是為了達到 邊緣增強；③非極大值抑制（匪S)。目的是保留梯度方向上的最大值。其貢獻在于給出了 一種尋找最優算子的思路，即確立了最優準則。
[0017] 步驟（3)、2D曲線在二維上是具有明顯視覺意義的線，其表征了體數據的拓撲和 幾何結構。為了更準確的表示體數據的結構信息將這些2D曲線反投回三維空間，找到每 條2D曲線對應3D候選曲線。然而2D曲線與3D曲線是一對多的關系，為了求解一條最有 3D對應曲線，定義了一個能量最小化函數，因此求解3D曲線問題轉化為一個全局最小化問 題。3D曲線上點有位置信息，不透明度信息，顏色信息，都是在最小化問題需要考慮的因子。 由于求解的最優3D曲線是曲線上的相連點之間能量最小，因此能量最小化問題可以轉化 為最短路徑問題，利用動態規劃算法可以求解該問題，并且借助GPU并行計算能力，大大提 高計算效率。在GPU上每個線程負責一條3D線的計算過程。
[0018] 2)曲線方向場生成
[0019] 步驟（4)、將2D曲線反投射回3D空間后，對3D空間進行采樣，求得的3D曲線在采 樣空間內會經過若干采樣點，利用卷積曲面的方法，對經過采樣點一定范圍的曲線進行積 分求得采樣點的曲線密度，以構建曲線密度場。曲線密度場中每個采樣點的密度值是衡量 其空間范圍內的曲線的多和少，由于很多采樣點周圍并沒有曲線經過，因此其曲線密度值 為零，對這些點則可以通過預計算剔除，只需計算剩余采樣點的曲線密度值即可。該過程是 在GPU上實現，在預計算階段，每個線程控制一條曲線，將其經過范圍的采樣點標記，沒被 標記的點則是被剔除的點。在采樣點曲線密度求解階段，每個線程控制一個剩余采樣點的 計算，大大提高了計算速度。
[0020] 步驟（5)、構建完曲線密度場后，用Sobel算子求得采樣點的梯度，由于步驟（4)已 經剔除了值為零的采樣點，所以此過程只需計算非零采樣點的梯度值。梯度表征了曲線密 度場的變化率最大的方向，因此采樣點的梯度大致垂直曲線方向，然后在利用最小二乘方 法可以在梯度基礎上求得采樣點的曲線方向，該方向與3D曲線的走勢基本一致。該過程是 在GPU上實現，每個線程控制一個非零采樣點的計算過程，大大提高了計算速度。
[0021] 3)體數據模型骨架的生成
[0022] 步驟（6)、生成曲線方向場后，由于空間采樣點與體素在空間位置上有對應關系， 將非零采樣點的方向賦值到其對應的體素上，以具有方向的體素的位置和方向為依據構建 包圍盒，包圍盒的構建過程是以連通性為依據，在體素方向上上下一定距離L處構建兩個 個垂直于體素方向的平面，兩個平面之間所有聯通的體素歸為一組，然后將包圍盒內所有 體素求一個平均值并作歸中處理，得到表征節點。
[0023] 步驟（7)、以步驟（6)得到的表征節點為節點，以表征節點之間的距離為權，構建 無向完全圖G，并用Prim算法求解圖G的最小生成森林，檢測出森林中樹的葉節點和分叉 點，從而進一步生成Polyline曲線。對于得到的Polyline曲線，由于是把各個表征節點點 進行連接得到，這樣的曲線是不夠圓滑的，采用B樣條擬合方法，Polyline上的點作為控制 點，計算出更多的插值點，求得B樣條曲線，使得曲線更光滑，得到最終體數據模型骨架。
[0024] 本發明與現有技術相比的優點在于：本發明能夠快速地提取出三維體數據模型的 骨架并且能得到較好的效果。本發明主要有兩點貢獻：第一，給出了一種基于渲染求解曲線 方向場的方法，該方法是從RayCasting直接體繪制算法渲染的2D圖像出發，提取其輪廓和 部分關鍵特征線，進而投射回3D空間，并根據3D候選曲線利用卷積曲面方法構建曲線密度 場，在此基礎上進而生成了曲線方向場，所有過程均在GPU上并行進行，具有快速高效并較 準確的優點。第二，將曲線方向賦值到體素上，并有體素位置和方向構建包圍盒求得表征節 點，連接這些節點并利用B樣條進行平滑擬合得到體數據模型的骨架，該過程降低了計算 量使得骨架提取過程快速而準確。

【專利附圖】

【附圖說明】
[0025] 圖1為本發明整體流程圖；
[0026] 圖2為RayCasting直接體繪制原理圖；
[0027] 圖3為射線穿越體數據示意圖；
[0028] 圖4為體繪制渲染的2D圖像效果圖；
[0029] 圖5為提取的2D曲線效果圖；
[0030] 圖6為2D曲線反投射回3D空間示意圖；
[0031] 圖7為3D候選曲線效果圖；
[0032] 圖8為計算采樣點曲線密度場示意圖；
[0033] 圖9為采樣空間曲線密度場效果圖；
[0034] 圖10為采樣點梯度效果圖；
[0035] 圖11為采樣點曲線方向效果圖；
[0036] 圖12為體素方向效果圖；
[0037] 圖13為表征節點效果圖；
[0038] 圖14為體數據模型骨架效果圖。

【具體實施方式】
[0039] 下面結合附圖以及本發明的【具體實施方式】進一步說明本發明。
[0040] 如圖1所示，給定體數據由體素組成，本發明針對三維體數據的骨架提取包括以 下步驟：
[0041] 步驟（1)如圖2(只設了 3個）在X、Y、Z3個坐標軸的正負方向（上下前后左右） 上設置六個相機，用RayCasting直接體繪制算法對三維體數據進行渲染，得到6幅2D圖 像；基于圖像空間的算法，對圖像的每一個像素，沿固定方向（通常是視線方向）發射一條 光線，光線穿越體數據的整個圖像序列，并在這個過程中，對圖像序列進行采樣獲取顏色信 息，同時依據光線吸收模型將顏色值進行累加，直至光線穿越整個圖像序列，最后得到的顏 色值就是渲染圖像的顏色。采用了從前面到背面進行采樣合成，其中顏色和不透明度的累 加公式如公式（1):
[0042] Cj = (l-α 卜廣)。，。-/ (1)
[0043] α J = (1-α 卜廣）α ηΛ
[0044] 其中，Q和α i分別是在體紋理上采樣所得到的顏色值和不透明度，其實也就是體 素中蘊含的數據；C/和α /表示累加的顏色值和不透明度。i表示當前步驟，i-Ι表示前 一步驟。
[0045] 如圖3所示，假定光線從f點投射到立方體中，并從1點投出，在立方體中穿越的 距離為m。當光線從f點投射到立方體中，穿越距離為n (n〈m)時進行采樣，則存在公式（2):
[0046] t = tstart+d* Δ (2)
[0047] 其中tstart表示立方體表面被投射點f的體紋理坐標；d表示投射方向；Λ表示采 樣間隔，隨著η的增加而遞增；t為求得的采樣紋理坐標。通過求得的采樣紋理坐標就可以 在體紋理上查詢體素數據。直到n>m，或者透明度累加超過1，一條射線的采樣過程才結束。
[0048] 該過程可以在GPU上實現，利用CUDA (計算統一設備架構）的硬件加速，大大提高 渲染速度。得到6個2D渲染結果，如圖4。
[0049] 步驟（2)、對步驟⑴得到的6幅2D圖像，使用輪廓檢測算法提取2D圖像的輪廓， 直接提取的輪廓需要收縮以獲得有效輪廓線，另外為表示更多細節需要使用Canny算子對 2D圖像進行部分關鍵特征線的提取。Canny算子是一階算子，其邊緣檢測可以分為三個步 驟：①高斯平滑函數。目的是為了平滑以消除噪聲；②一階差分卷積模板。目的是為了達到 邊緣增強；③非極大值抑制（匪S)。將提取的輪廓線與用Canny提取出特征線融合得到效 果圖圖5,圖5中同一顏色的線代表為一 2D曲線。
[0050] 步驟（3)、2D曲線在二維上是具有明顯視覺意義的線，其表征了體數據的拓撲和 幾何結構。為了更準確的表示體數據的結構信息，將這些2D曲線反投回三維空間，如圖6其 中由空心的點連接成的曲線是一條2D曲線，需要找到2D曲線對應3D候選曲線，在圖6中由 實線連接黑色實心點而成的曲線是得到的3D曲線。由圖6可知一個2D點反投射回3D空間 對應多個3D點，即圖6中一個空心點對應多個實心點。定義了一個能量最小化函數（3):

【權利要求】
1. 一種基于渲染的快速體數據骨架提取方法，其特征在于包括如下步驟： (1) 在X、Y、Z，3個坐標軸的正負方向即上下前后左右上設置六個視點，用RayCasting 直接體繪制算法對三維體數據進行渲染，得到6個2D圖像； (2) 利用輪廓線檢測方法對得到的6個2D圖像進行輪廓線的提取時，直接提取的輪廓 需要收縮以獲得有效輪廓線，并用Canny算子提取所有2D圖像的部分關鍵特征線，將輪廓 線與部分關鍵特征線融合以獲得有效2D曲線，這些2D曲線體現了體數據的幾何及拓撲結 構；所述部分關鍵特征線即關節部位的特征線； (3) 對步驟（2)提取的2D曲線反投射回三維，由于每條2D曲線與3D曲線是一對多的 關系，故定義一個能量最小化函數，利用動態規劃算法求解最優解，得到2D曲線所對應的 最優3D曲線；該過程是在GPU上實現，每個線程控制一條曲線的求解過程，大大提高了求解 速度；該過程是在GPU上實現，每個線程控制一條曲線的求解過程，大大提高了求解速度； (4) 對空間進行采樣，求得的3D曲線在采樣空間內會經過若干采樣點，采用卷積曲面 的方法，對經過采樣點一定范圍的曲線進行積分求得采樣點的曲線密度，以構建曲線密度 場；由于很多采樣點周圍并沒有曲線經過，因此其曲線密度值為零，對這些點則通過預計算 剔除，只需計算剩余采樣點的曲線密度值即可；該過程是在GPU上實現，在預計算階段，每 個線程控制一條曲線，將其經過范圍的采樣點標記，沒被標記的點則是被剔除的點；在采樣 點曲線密度求解階段，每個線程控制一個剩余采樣點的計算，大大提高了計算速度； (5) 獲得曲線密度場后，用Sobel算子求得采樣點的梯度，由于步驟（4)已經剔除了值 為零的采樣點，所以此過程只需計算非零采樣點的梯度值；梯度表征了場的變化率最大的 方向，因此采樣點的梯度大致垂直曲線方向，利用最小二乘方法可以在梯度基礎上求得采 樣點的曲線方向，該過程是在GPU上實現，每個線程控制一個非零采樣點的計算過程，大大 提高了計算速度； (6) 由于3D曲線表征了體數據的幾何和拓撲結構，求得采樣點的曲線方向后，將采樣 點的曲線方向賦值給其所對應的體素，以體素的位置和方向為依據構建包圍盒，包圍盒的 構建過程是以連通性為依據，將包圍盒內所有體素求平均值，得到表征節點； (7) 用Prim算法將所用表征節點用最小生成樹連接成Polyline曲線，最后用B樣條將 求得的Polyline曲線平滑擬合，得到最終體數據骨架。
2. 根據權利要求1所述的基于渲染的快速體數據骨架提取方法，其特征在于：所述步 驟（1)包括：從X，Y，Z正負軸方向設置6個相機，用RayCasting直接體繪制算法對三維體 數據進行渲染，得到6個2D圖像，該過程是在GPU上實現，利用計算統一設備架構CUDA的 硬件加速，大大提高渲染速度，并利用PBO, FB0技術渲染到紋理，并保存到2D渲染結果。
3. 根據權利要求1所述的基于渲染的快速體數據骨架提取方法，其特征在于：所述步 驟⑵具體實現為：以這些表征節點為節點，以表征節點之間的距離為權，構建無向完全圖 G，并用Prim算法求解圖G的最小生成森林，檢測出森林中樹的葉節點和分叉點，從而進一 步生成Polyline曲線。采用B樣條擬合方法，Polyline上的點作為控制點，計算出更多的 插值點，求得B樣條曲線，得到更加的光滑體數據骨架。
【文檔編號】G06T15/00GK104156997SQ201410363994
【公開日】2014年11月19日 申請日期:2014年7月28日 優先權日:2014年7月28日
【發明者】王莉莉, 戚慶林, 侯飛, 李帥, 郝愛民 申請人:北京航空航天大學