基於最大公共子圖與包圍盒的斷骨斷面及斷骨模型配準方法與流程

2023-04-23 00:50:26

本發明涉及一種基於最大公共子圖與包圍盒的斷骨斷面及斷骨模型配準方法。涉及專利分類號g06計算；推算；計數g06t一般的圖像數據處理或產生g06t7/00圖像分析，例如從位像到非位像。

背景技術：

目前骨折手術普遍採用人工復位和傷肢內固定相結合的方法，這種方法存在的問題是創傷大、出血多、容易引發神經血管損傷等併發症。因此，可以利用計算機對斷骨模型進行虛擬拼接，從而在術前得到鋼板的各種幾何參數。然而目前的斷骨虛擬拼接方法存在很多問題，比如斷面分割不準確、需要手動操作、預配準不精確等等，這些問題都很大程度地影響了斷骨模型配準的效果。

在現有技術中，存在通過根據斷面法向量與主軸線夾角和法向量的突變兩個方面來進行斷面分割的手段，也存在直接根據斷骨的軸線進行預配準的方案。但是，根據斷面法向量與主軸線夾角和法向量的突變兩個方面來進行斷面分割，需要人為的在模型上進行點選，不能實現斷面的自動提取。

另外，單純將斷骨的主軸線進行對齊並不能滿足預配準的要求，可能會出現斷骨的截面方向相背離的情況，而且截面形狀也無法做到大致吻合，很容易造成精配準的不準確。

技術實現要素：

本發明要解決的技術問題是：一種基於最大公共子圖與包圍盒的斷骨斷面及斷骨模型配準方法，包括如下步驟：

斷骨軸線提取步驟：

根據經由ct圖像序列重建生成的斷骨模型，生成該模型初始樣本矩陣z；對該樣本矩陣進行中心化，得到矩陣x；計算該矩陣x的協方差矩陣c和協方差矩陣c的特徵值，選取最大特徵值對應的特徵向量作為斷骨軸線；

斷骨高斯映射步驟：

計算所述斷骨模型上每一個三角面片的法向量，對全部的法向量進行單位化，得到每一個三角面片的單位法向量n；將得到的所有單位法向量n的起點平移到坐標原點，向量的終點則會落到以坐標原點o為球心，半徑為1的單位球面s上，斷骨模型的高斯映射完成；

斷面點集提取步驟：

通過剔除其法向量與所述斷骨軸線呈一定角度內的三角面片，將所述的斷骨模型分割成兩個斷骨模型；每個模型包含該斷骨的斷面模型；

分別對兩個斷骨進行三角面片的剔除，最終得到四個子模型，其中2個為斷面模型；

提取4個所述子模型包含節點數的斷面點集，通過使用最大公共子圖算法兩兩比較所述的斷面電機，選擇2個包含節點數最大的一對子模型作為所述的兩個斷面模型；

斷骨預配準步驟：

建立所述的兩個斷面模型的包圍盒，將兩包圍盒在空間上對齊完成預配準；

斷骨精確配準步驟：

通過icp算法對兩斷骨模型進行精確配準，完成兩斷骨的接合。

作為優選的實施方式，所述的斷骨軸線提取步驟具體包括如下步驟：

a.對斷骨模型中的n個頂點，將所述的n個點坐標賦值給一個3×n的矩陣z，生成該模型的初始樣本矩陣為三維列向量；

b.對所述的初始樣本矩陣進行中心化，得到中心化後的矩陣x；

其中：計算得到的三維向量則可以視為斷骨模型在空間上的中心；

c.計算所述矩陣的協方差陣列c，轉換過程如下公式所示：

計算所述的協方差陣列c為3×3矩陣；

d.分別求出所述協方差陣列c的特徵值λ1，λ2，λ3和特徵值對應的特徵向量

比較所述的特徵值λ1，λ2，λ3，選擇最大的特徵值對應的特徵向量，作為該斷骨模型的軸線向量。

作為優選的實施方式，所述的斷骨模型的高斯映射步驟具體包括如下步驟：

a.計算所述三維模型上每個三角面片的法向量；

b.對得到的所述每個三角面片的法向量進行單位化，得到每個所述三角面片的單位法向量n；

c.將得到的每個單位方向量的起點平移至坐標原點，則每個所述的單位法向量的重點最終落在以坐標原點o為球心，半徑為1的單位球面s上，完成高斯映射。

作為優選的實施方式，所述的斷面點集提取步驟中剔除的過程如下：

計算每個所述三角面片法向量與所述斷骨模型的軸線方向向量的夾角；剔除法向量所成夾角在0.4π-0.6π之間的三角面片，將剩餘的法向量與軸線方向向量夾角範圍在[0,0.4π]∪[0.6π,π]內的三角面片保留，將一個斷骨模型分割成2個子模型，其中一個所述的子模型為斷面模型；

對所述的2個斷骨模型中將夾角在0.4π-0.6π之間的三角面片剔除，完成斷骨模型側面的三角面片剔除，得到四個子模型，其中兩個是斷面模型；

更進一步的，得到兩個斷面模型的過程如下：

使用最大公共子圖算法判斷斷面模型，分別對一個斷骨模型的兩個子模型與另一斷骨的兩個子模型使用最大公共子圖算法，最大公共子圖所包含節點數最大的一對子模型，即為2個斷面模型。

更進一步的，所述的精確配準斷骨精確配準步驟過程如下：

將兩斷面模型分點集分別記為p與q；在目標點集p中選取點集合pik∈p，計算源點集合q中相應的點滿足

計算旋轉矩陣rk和平移矢量tk，滿足計算pk+1＝{pik+1|pik+1＝rkpik+tk,pik∈p}和其中，k代表著第k次迭代，pik和分別代表目標點集和源點集中的點，min代表設定好的最小距離的閾值，dk+1代表著第k+1次迭代下的平均最小距離；

若dk+1大於等於事先設定的閾值t，則重新在目標點集p中選取點集合步驟，重新開始算法，當dk+1<t或循環次數大於事先設定好的循環次數的閾值時，跳出循環。

作為優選的實施方式，斷骨精確配準步驟之後還具有虛擬鋼板預彎步驟：

在拼接好的斷骨斷面附近進行點選，確定鋼板模型的形狀以及尺寸；

記錄下點選的三角平面值，選中範圍內的所有表面三角面片；

計算每個三角面片的法向量值；將每個平面根據其法向量方向進行一定程度的加厚，並填充其縫隙的部位；

所得到的加厚部分即為模擬的鋼板模型三維數據。

附圖說明

為了更清楚的說明本發明的實施例或現有技術的技術方案，下面將對實施例或現有技術描述中所需要使用的附圖做一簡單地介紹，顯而易見地，下面描述中的附圖僅僅是本發明的一些實施例，對於本領域普通技術人員來講，在不付出創造性勞動的前提下，還可以根據這些附圖獲得其他的附圖。

圖1為本發明的算法流程圖

圖2為本發明斷骨軸線提取後的效果示意圖

圖3為本發明斷骨模型高斯映射的結果示意圖

圖4為本發明斷骨模型分割成4個子模型的分割結果示意圖

圖5為本發明使用最大公共子圖方法提取的斷面模型示意圖

圖6為本發明基於包圍盒對齊的斷骨模型預配準結果示意圖

圖7為本發明斷骨模型精確配準結果示意圖

圖8為本發明斷骨表面選取的特徵點示意圖

圖9為本發明虛擬鋼板模型示意圖

具體實施方式

為使本發明的實施例的目的、技術方案和優點更加清楚，下面結合本發明實施例中的附圖，對本發明實施例中的技術方案進行清楚完整的描述：

如圖1-9所示：基於最大公共子圖與包圍盒的斷骨斷面及斷骨模型配準方法，主要包括如下步驟：

斷骨軸線提取

本方法利用pca算法進行斷骨的軸線提取。提取效果如圖2所示。

輸入初始樣本矩陣。根據輸入的模型，設模型一共包含n個點，將這n個點的坐標值賦值給一個3×n的矩陣，此矩陣則為初始樣本矩陣

初始樣本矩陣中心化。將初始樣本矩陣z中心化後，得到矩陣中心化的公式為其中，計算得到的三維向量則可以視為斷骨模型在空間上的中心。

計算矩陣x的協方差陣列c。c的公式為計算得到的c是一個3×3的矩陣。

求協方差矩陣c的特徵值和特徵向量。對於得到的協方差陣列c，求出其特徵值λ1，λ2，λ3與對應的特徵向量將三個特徵值做比較，得到最大的特徵值，其對應的特徵向量，就是該斷骨模型的軸線方向。

斷骨高斯映射

本方法將斷骨模型進行了高斯映射操作，映射結果如圖3所示。計算斷骨模型上的每一個三角面片的法向量。

將計算得到的三角面片的法向量單位化，單位化的過程可採用如下所述的方法：

設法向量為e＝(e1，e2，e3)，那個其單位法向量

得到每個三角面片的單位法向量n。將得到的所有單位法向量的起點平移到坐標原點，向量的終點則會落到以坐標原點o為球心，半徑為1的單位球面s上，斷骨模型的高斯映射完成。

通過在高斯映射中將每個面片的法向量計算出來，通過對高斯映射的顯示，可以直觀的觀察到法向量與模型軸線的關係，為下一步根據法向量與軸線夾角的關係剔除三角面片做鋪墊。

斷骨斷面數據點集的提取

斷骨模型分割。根據每個三角面片的法向量，計算其與模型軸線方向向量的夾角，將與軸線方向向量夾角在(0.4π,0.6π)範圍內的三角面片剔除，則留下的三角面片法向量與軸線方向向量的夾角範圍是[0,0.4π]∪[0.6π,π]，這樣，就把一個斷骨模型分割成了兩個子模型，其中包含該斷骨的斷面模型。分別對兩個斷骨進行三角面片的剔除，則可以獲得四個子模型，其中有兩個是斷面模型。模型分割結果如圖4所示。

基於最大公共子圖的斷面數據點集提取。本方法提出了一種基於最大公共子圖算法的方法來提取斷面點集。提取結果如圖5所示。本方法的大致實現步驟是：

構建圖結構。將模型的一個三角面片看作是圖結構中的一個頂點，如果兩個三角面片在模型中是邊相鄰的關係，那麼代表兩個三角面片的頂點之間則存在一條邊將這兩個頂點連接起來。用該頂點所表示的三角面片的法向量與其邊相鄰的三個三角面片的法向量的夾角之和來表示每個頂點的值。對於每條邊的值，只要兩個頂點之間有邊的關係，則將該邊的值賦值為1。

使用最大公共子圖算法判斷斷面模型。分別對每個子模型與另一個模型的兩個子模型使用最大公共子圖算法，最大公共子圖所包含的節點數最大的一對子模型，即為兩個斷面模型。

最大公共子圖算法流程如下：

(1)對於已經建立好的兩個圖g1、g2，構建關聯圖g；

(2)根據圖g的結構，分別對remaining\connected和remaining∩connected兩個節點集進行染色，得到colourclasses；

(3)對colourclasses中每個顏色類中的每一個節點v進行判斷，如果當前的結果集solution的節點個數與colourclasses的類的個數小於最優解incumbent的節點個數或者當前節點不在connected點集裡並且solution集合不為空集，則刪除當前顏色類，從下一個顏色類進行節點遍歷；

(4)否則，將當前節點v添加進solution集合裡，如果新的solution的節點個數大於目前最優解incumbent的節點個數，則更新incumbent集合，並更新connected和remaining集合；

(5)如果當前remaining集合不為空，則向下遞歸尋找最大公共子圖，直到remaining集合為空集。

對於以上流程，其中，solution集合中包含當前已經尋找到的公共子圖的節點，remaining集合中包含與solution中的節點相鄰的和還沒有被接受或者拒絕的節點，connected集合中包含在g1中與已經被接受的頂點相鄰的頂點所匹配的節點對。

基於包圍盒對齊的斷骨預配準

斷骨模型的預配準，是通過將兩個斷面模型的包圍盒在空間上對齊的方式來實現的。預配準結果如圖6所示。

將包圍盒的短軸方向固定在y軸上，長軸方向固定在x軸上，中軸方向固定在z軸方向。如果發現某一個包圍盒的長軸方向與中軸方向分別在z軸和x軸上，則需要對該包圍盒以及對應的斷骨模型進行空間變換。先是將包圍盒以及斷骨模型繞y軸旋轉90度，再沿著x軸方向平移包圍盒長軸的長度，這樣就可以實現使包圍盒三個軸的方向分別於三個坐標軸對應。

經過第一步變換之後，斷骨模型會有兩種狀態，一種是整個斷骨模型完全在xoz平面的一側，另一種是斷骨模型被xoz平面分為兩部分。對於第二種情況，我們還需要對模型位置進行處理，將其平移到xoz平面的一側。對於平移的方向與距離，如果斷骨模型的中心點的y坐標大於0，則將模型沿著y軸正方向平移包圍盒短軸長度的距離，如果斷骨模型的中心點的y坐標小於0，則將模型沿y軸負方向平移包圍盒短軸長度的距離。

兩個斷骨模型被變換到了xoz平面的同側或異側。對於在xoz平面同側的兩個斷骨模型，需要將其變換到分別在xoz平面的兩側。變換的過程為，將一個斷骨模型繞著x軸或者z軸旋轉180度，再沿著z軸方向或者x軸方向平移中軸或者長軸長度的距離。

將兩個斷骨模型的中心點坐標與包圍盒的中心點坐標進行比較，如果包圍盒中心點的x坐標在兩個斷骨模型中心點坐標的中間，則兩個斷骨已經預配準完成；如果包圍盒中心點的x坐標比兩個斷骨模型中心點的x坐標都大或者都小，則說明發生了180度翻轉的情況。此時，只需要將一個模型繞著y軸旋轉180度，並在x軸方向平移包圍盒長軸長度，在z軸方向平移包圍盒中軸長度即可。

斷骨精配準

預配準結束後，利用icp算法對兩組斷骨模型點集進行精配準操作，進一步提高其準確性。精配準結果如圖7所示。

將兩部分點集分別記為p與q。

在目標點集p中選取點集合pik∈p，計算源點集合q中相應的點滿足

計算旋轉矩陣rk和平移矢量tk，滿足

計算pk+1＝{pik+1|pik+1＝rkpik+tk,pik∈p}和

若是dk+1大於等於事先設定的閾值t，則跳轉到步驟2，當dk+1<t或循環次數大於事先設定好的循環次數的閾值時，跳出循環。

虛擬鋼板預彎

在本步驟，主要利用對斷骨斷裂部位曲面數據的相應操作來模擬出鋼板的大概形態。

在拼接好的斷骨斷面附近進行點選(如圖8)，確定鋼板模型的大概形狀以及大小。

記錄下點選的三角平面值，選中範圍內的所有表面三角面片。計算每個三角面片的法向量值，並記錄下來。將每個平面根據其法向量方向進行一定程度的加厚，並填充其縫隙的部位。所得到的加厚部分即為模擬的鋼板模型三維數據(如圖9)，可將其作為結果導出輸出。

以上所述，僅為本發明較佳的具體實施方式，但本發明的保護範圍並不局限於此，任何熟悉本技術領域的技術人員在本發明揭露的技術範圍內，根據本發明的技術方案及其發明構思加以等同替換或改變，都應涵蓋在本發明的保護範圍之內。