本發(fā)明是關于圖像處理技術領域，特別涉及一種肋骨可視化的肋骨骨折輔助診斷方法。

背景技術：

計算機斷層掃描(ct)是診斷肋骨骨折的主要方法，但肋骨特殊的組織結構導致肋骨遮擋因而微小骨折及骨裂容易漏診，精確的診斷需要經驗的積累，常規(guī)的二維方法需要手動追蹤多個切片，傳統(tǒng)的三維可視化方法在診斷時需要對圖像進行旋轉觀察從而避免肋骨遮擋的影響。常用重建方法有多層面重組法、最大密度投影法及表面陰影顯示法。多層面重組法需要做一系列圖像且非目標器官存在形變；最大密度投影法低密度組織器官易被遺漏且噪聲影響大；表面陰影顯示法容易產生偽影且看不出密度信息。

技術實現(xiàn)要素：

本發(fā)明的主要目的在于克服現(xiàn)有技術中的不足，提供一種通過把肋骨投影到二維平面上，消除在臨床診斷時肋骨遮擋導致微小骨折及骨裂的骨折輔助診斷方法。為解決上述技術問題，本發(fā)明的解決方案是：

提供一種肋骨可視化的肋骨骨折輔助診斷方法，包括下述過程：

一、獲取圖片；

二、提取肋骨；

三、人工修正，去掉脊柱、胸骨；

四、展開肋骨；

所述過程一是指：取得原始的人體胸部ct掃描圖像；

所述過程二是指：將人體胸部ct掃描圖像轉化為256階灰度圖像，然后用增強濾波器對肋骨進行增強處理，調整像素值大小取一個范圍提取出整個胸部骨骼，整個胸部骨骼包括肋骨、脊柱及胸骨，按照肋骨空間位置大小從上向下標記左右肋骨分支；

所述過程三是指：創(chuàng)建三維窗口，用戶通過鼠標右鍵點擊劃圈標示出不屬于肋骨的部分，計算機根據用戶的標示將圈內部分(脊柱、胸骨)去除；

所述過程四是指：泊松曲面重建肋骨輪廓，拉普拉斯網格收縮提取肋骨骨架；應用最近點迭代矯正骨架，對準骨架位置；應用曲面重建算法(curvedplanarreformation，cpr)對每一只肋骨進行展開，得到最終的肋骨可視化展開結果，用于輔助醫(yī)生進行肋骨病變診斷。

在本發(fā)明中，所述過程二具體包括下述步驟：

步驟a：把人體胸部ct掃描圖像轉化為256色灰度圖像

步驟b：然后使用一個增強濾波算子對肋骨區(qū)域進行強化，從圖像中逐點提取平滑的黑塞矩陣：

其中，所述i是指每點的二階導數矩陣；所述iij是指二階導數，i、j分別取x、y或者z；

步驟c：黑塞矩陣i通過高斯核的二階偏導數與原圖像進行卷積求得：

其中，gσ(v)＝exp(-v^tv/σ²)，v∈r³為高斯核；所述δ為拉普拉斯算符；所述為步驟a中的256色灰度圖像；所述σ是指當前高斯核的標準差；所述t是指向量的轉置；所述exp是指自然指數函數；所述*指的是函數間的卷積；所述r³為三維歐氏空間；所述v為歐氏空間中的一個三維向量；

步驟d：求出黑塞矩陣i的特征值λ1，λ2，λ3，且|λ1|＜|λ2|＜|λ3|；

計算每點的增強值：

其中，所述v0(i)用于表示像素值類似血管(vesselness)的程度；所述α、β、c都分別大于0(均為可調整經驗參數，可根據實際灰度圖像自行調節(jié))；

步驟e：通過調整像素值范圍提取出胸部骨骼，包含肋骨、脊柱及胸骨；

步驟f：按照肋骨空間位置大小從上到下對左右肋骨進行標記。

在本發(fā)明中，所述過程四中，泊松曲面重建肋骨輪廓具體為：

輸入經過過程三人工修正后提取的肋骨分割結果的表面點集，記為點集合s，該點集合滿足任意元素s∈s，且任意元素s都具有兩個已知參數：元素s的三維空間中對應的位置s.p、元素s的內法線方向同時，記vec為三維歐式空間上的向量場，該向量場滿足對所有的s∈s，向量場的值在點s.p上的值為而向量場在其余點的值為0：

泊松重建是尋找三維空間上的一個標量函數χ，使得式子最?。黄渲?，符號||·||代表歐氏距離范數；

泊松重建轉而求解將作為上述χ的近似解；上式中，△代表拉普拉斯算符，為三維空間上的一個標量函數；

最后，泊松重建通過求取同屬一個水平集的點集：獲得表面重建結果的頂點集合vp，并對vp構造相應八叉樹，即獲得對應表面重建結果gp，輸出的泊松重建的結果記為：gp＝(vp，ep)；

上式中，vp代表待求水平集上的點集，q代表三維歐氏空間中的點，其中|·|代表集合的元素個數；gp代表泊松重建表面網格結果，其中ep是表面網格中邊的集合，ep中任意元素是vp中兩點間的連線。

在本發(fā)明中，所述過程四中，拉普拉斯網格收縮提取肋骨骨架是指：在提取骨架過程中保持了細節(jié)形變，通過估計肋骨上點云鄰域建立散亂點之間的連接關系，然后構建拉普拉斯加權矩陣，將所有點作為全局的位置約束，進行點云收縮，對收縮后的點云進行聚類，連接成線得到肋骨骨架；具體為：

對于任意一個表面重建結果g＝(v，e)；其中，v為該表面的頂點集合(泊松重建所提取的等值面中的點)；e是該表面重建結果中邊的集合(表面網格中邊的集合)，是v中點的連線；假設g中的頂點數，即集合v的元素個數為n，記v＝{v1，...，vn}，vi∈r³；定義并計算該表面重建結果的拉普拉斯矩陣為l，對于該表面重建結果g(v，e)，記αij和βij為邊(i，j)∈e的對角，其中(i，j)∈e代表從vi到vj的連線包含在e中，i，j分別為v中點的索引；定義角度權重：

wi，j＝cotαi，j+cotβi，j

通過如下方式計算局部拉普拉斯矩陣l＝{li，j}n×n

上式中，(i，j)中的i和j分別指代當前邊的兩個頂點的索引；k指代在邊(i，k)中不同于當前索引i的那個頂點的索引；li，j是指拉普拉斯矩陣第i行第j列的元素；

拉普拉斯網格收縮提取肋骨骨架利用一個迭代過程實現(xiàn)，記上一步驟中得到的肋骨泊松表面重建結果為gp＝g¹＝(v¹，e)，其中g¹，v¹和e分別代表算法初始時的表面重建結果、與之對應的初始頂點集合以及邊的集合，v¹＝vp，e＝ep；設初始權重矩陣和

首先計算對應的拉普拉斯矩陣l¹，在第t次迭代過程中，通過已知的第t步的結果其變量分別為第t步的表面重建結果，權重，以及拉普拉斯矩陣，算法重復以下步驟，并更新結果到第t+1次迭代中，直到當前迭代生成的表面重建結果體積為0，迭代過程具體為：

1)更新第t+1步的表面重建結果g^t+1(v^t+1，e)，其中v^t+1是第t+1次更新后的頂點集合，通過求解最小二乘：[wl^tl^twh^t]^tv^t+1＝[0wh^tv^t]^t；

其中，大寫的t代表向量的轉置，v^t為第t次迭代生成的頂點集合，wl^t以及wh^t為第t次迭代的權重矩陣；

2)更新第t+1步的權重矩陣為：

其中，和分別是頂點i當前鄰域的面積和最初鄰域的面積，sl值可調(一般取2.0)；

3)通過表面重建結果g^t+1(v^t+1，e)，計算第t+1步的拉普拉斯矩陣l^t+1；

迭代過程結束后，最后生成的表面重建結果為g^∞(v^∞，ep)，并且獲得該最后收斂的表面重建結果的骨架v^∞，v^∞為該表面的頂點集合。

在本發(fā)明中，所述過程四中，應用最近點迭代矯正骨架具體是指：對給定兩組不同的三維數據點集，其中一個為源數據點集(待配準點集)p＝{p1，..，pn}，另一個為目標點集q＝{q1，...，qm}，其中m≥n，icp試圖找出一種剛性變換，將點集p與目標q進行匹配，即icp通過迭代的方法來尋找這樣的最佳剛性匹配；

在第k次迭代過程中，icp算法找到與源數據點集p中的每個點對應的目標點集q中的點，并且計算最佳的剛性變換，使得下述距離最?。?/p>

其中，分別為第k步迭代的點集p，q中的對應點，t＝(r；t)，r∈so(3)，t∈r³為剛性變換，γ代表所有剛性變換的集合；其中so(3)代表三階正交群，||·||代表歐氏度量，符號代表函數復合；

具體來說icp主要包括下述步驟：

(1)對每一個pi^k，計算目標點集q中對應點使最?。?/p>

(2)計算旋轉矩陣r^k與平移向量t^k，極小化距離：

(3)更新第k+1步源點云中點的位置：

(4)計算l²距離誤差：

(5)如果d^k+1不小于給定的閾值δ，返回到(1)執(zhí)行，直到d^k+1＜δ或迭代次數大于預設的最大迭代次數為止；

最終，得到的變換為：其中ti為第i次迭代中求得的剛性變換；

因源數據點集p為肋骨骨架內的骨骼點v^∞，而點集q為與之對應的肋骨二值分割結果(看作點集)，通過上述應用最近點迭代矯正骨架的迭代算法，即實現(xiàn)將肋骨骨架與分割結果對齊。

在本發(fā)明中，所述過程四中，cpr方法展開肋骨具體包括下述步驟：

(1)應用得到的肋骨骨架(即通過應用最近點迭代矯正后的骨架)，均勻選取上面的3d采樣點；

(2)在每個采樣點上計算基于相鄰點和前一點的正交觀測方向，在每個點附近采樣一個平面；

(3)在每個點上重復上述過程，取一個切片，重復第二步；

(4)把采樣的切片相互堆疊在一起形成一個直卷；

(5)重復每根肋骨從而把肋骨展開。

與現(xiàn)有技術相比，本發(fā)明的有益效果是：

本發(fā)明通過把肋骨展開投影到平而上，消除了肋骨彎曲遮擋對診斷的影響，相對與傳統(tǒng)臨床診斷方法，在平面上觀察提高率診斷的直觀性，本發(fā)明可以作為診斷微小骨折骨裂的輔助方法。

附圖說明

圖1為肋骨二值圖。

圖2為修正后去掉脊柱和胸骨。

圖3為肋骨xz方向展開圖。

圖4為肋骨yz方向展開圖。

具體實施方式

下面結合附圖與具體實施方式對本發(fā)明作進一步詳細描述：

一種新的肋骨可視化肋骨骨折輔助診斷方法，包括下述過程：

一、獲取圖片；

二、提取肋骨；

三、人工修正，去掉脊柱、胸骨；

四、展開肋骨。

所述過程一具體包括：通過和醫(yī)院聯(lián)系得到人體胸部ct掃描圖片。

所述過程二具體包括：

對腹部ct圖像進行處理，是基于骨骼的密度比較大，在ct掃描中ct值比較大，表現(xiàn)為骨骼與其他組織器官相比亮度高，因而有明顯區(qū)別，通過選取適當的窗寬窗位使用多尺度的增強濾波算子提取出肋骨，在保留骨骼區(qū)域的同時避免了其他組織和器官(如血管、鈣化點、腎臟)的干擾。具體步驟是：

步驟a：把圖像轉化為256色灰度圖像

步驟b：然后使用一個增強濾波算子對目標區(qū)域進項強化，從圖像中逐點提取平滑的黑塞矩陣

其中，所述i是指每點的二階導數矩陣；所述iij是指二階導數，i、j分別取x、y或者z；

步驟c：黑塞矩陣i通過高斯核的二階偏導數與原圖像進行卷積求得：

其中，gσ(v)＝exp(-v^tv/σ²)，v∈r³為高斯核；所述△為拉普拉斯算符；所述為步驟a中的256色灰度圖像；所述σ是指當前高斯核的標準差；所述t是指向量的轉置；所述exp是指自然指數函數；所述*指的是函數間的卷積，所述r³為三維歐氏空間，所述v為歐氏空間中的一個三維向量；

步驟d：求出黑塞矩陣i的特征值λ1，λ2，λ3，且|λ1|＜|λ2|＜|λ3|；

計算每點的增強值：

其中，所述v0(i)用于表示像素值類似血管(vesselness)的程度；所述α、β、c都分別大于0(均為可調整經驗參數，可根據實際灰度圖像自行調節(jié))；

步驟e：通過調整像素值范圍提取出胸部骨骼，包含胸骨、脊柱；

步驟f：按照肋骨空間位置大小從上到下對左右肋骨進行標記。

所述過程三具體包括：

創(chuàng)建三維窗口，用戶通過鼠標右鍵點擊劃圈標示出不屬于肋骨的部分，計算機根據用戶的標示將圈內部分(脊柱、胸骨)去除。

所述過程四具體包括：

泊松曲面重建肋骨輪廓是通過估計模型的指示函數和提取等值面，對表面重建一個無縫的三角逼近；向量場v為空間m(在模型表面具有高分辨率，遠離表面地方具有低分辨率)中的函數線性總和，構造并解泊松方程，提取求取的指示函數的等值面(它是一種全局的方法，把有法向點集的表面重建轉化為一個空間的泊松問題，從而把表面重建問題表示為泊松方程的解)。

輸入經過過程三人工修正后提取的肋骨分割結果的表面點集，記為點集合s，該點集合滿足任意元素s∈s都具有兩個已知參數：該元素s的三維空間中對應的位置s.p，以及該元素的內法線方向同時，記vec為三維歐式空間上的向量場，該向量場滿足對所有的s∈s，向量場的值在點s.p上的值為而向量場在其余點的值為0。泊松重建希望尋找三維空間上的一個標量函數χ，使得式了最小。其中符號||·||代表歐氏距離范數。泊松重建轉而求解然后將作為上述χ的近似解。上式中△代表拉普拉斯算符，為三維空間上的一個標量函數。最后，泊松重建通過求取同屬一個水平集的點集：獲得表面重建結果的頂點集合vp，并對vp構造相應八叉樹并獲得對應表面重建結果gp。上式中vp代表待求水平集上的點集，q代表三維歐氏空間中的點，其中|·|代表集合的元素個數。gp代表泊松重建表面網格結果。記gp＝(vp，ep)，其中ep是表面網格中邊的集合，ep中任意元素vp中兩點間的連線，下同。

在提取骨架過程中保持了細節(jié)形變，通過估計肋骨上點云鄰域建立散亂點之間的連接關系，然后構建拉普拉斯加權矩陣，將所有點作為全局的位置約束，進行點云收縮，對收縮后的點云進行聚類，連接成線得到肋骨骨架。

對任意一個表面重建結果g＝(v，e)，其中v為該表面的頂點集合。而e是該表面重建結果中邊的集合，為v中點的連線。假設g中的頂點數，即集合v的元素個數為n。在此情況下，記v＝{v1，...，vn}，vi∈r³。在這樣的定義下，我們定義并計算該表面重建結果的拉普拉斯矩陣為l。對于該表面g(v，e)，記αij和βij為邊(i，j)∈e的對角，其中(i，j)∈e代表從vi到vj的連線包含在e中，i，j分別為v中點的索引。定義角度權重：

wi，j＝cotαi，j+cotβi，j

我們通過如下方式計算局部拉普拉斯矩陣l＝{li，j}n×n

上式中，(i，j)中的i和j分別指代當前邊的兩個頂點的索引。k指代在邊(i，k)中不同于當前索引i的那個頂點的索引。li，j是指拉普拉斯矩陣第i行第j列的元素。

拉普拉斯網格收縮提取肋骨骨架利用一個迭代過程實現(xiàn)，我們記上一步驟中得到的肋骨泊松表面重建結果為gp＝g¹＝(v¹，e)，其中g¹，v¹和e分別代表算法初始時的表面重建結果、與之對應的初始頂點集合以及邊的集合，v¹＝vp，e＝ep。設初始權重矩陣和算法首先計算對應的拉普拉斯矩陣l¹。在第t次迭代過程中，通過已知的第t步的結果其變量分別為第t步的表面重建結果，權重，以及拉普拉斯矩陣，算法重復以下步驟，并更新結果到第t+1次迭代中，直到當前迭代生成的表面重建結果體積為0：

1)更新第t+1步的表面重建結果g^t+1(v^t+1，e)，其中v^t+1是第t+1次更新后的頂點集合，通過求解最小二乘：[wl^tl^twh^t]^tv^t+1＝[0wh^tv^t]^t。其中大寫的t代表向量的轉置，v^t為第t次迭代生成的頂點集合，wl^t以及wh^t為第t次迭代的權重矩陣。

2)更新第t+1步的權重矩陣為：其中，和分別是頂點i當前鄰域的面積和最初鄰域的面積，sl值可調，一般取2.0。

3)通過表面重建結果g^t+1(v^t+1，e)計算第t+1步的圖像拉普拉斯矩陣l^t+1。

記最后生成的表面重建結果為g^∞(v^∞，ep)，則該表面的骨架即為v^∞，其中v^∞為該表面的頂點集合。

在本專利中，算法將第三步生成的表面重建結果gp＝(vp，ep)作為初始值即g¹＝(v¹，e)帶入上述算法，并且獲得該表面重建結果的骨架。這個骨架是一系列點的集合，該集合即為最后收斂的表面重建結果g^∞(v^∞，ep)中的v^∞。

應用最近點迭代矯正骨架具體是指：對給定兩組不同的三維數據點集，其中一個為源數據點集(待配準點集)p＝{p1，...，pn}，另一個為目標點集q＝{q1，...，qm}，其中m≥n，icp試圖找出一種剛性變換，將點集p與目標q進行匹配，即icp通過迭代的方法來尋找這樣的最佳剛性匹配)。

在第k次迭代過程中，icp算法找到與源點集p中的每個點對應的目標點集q中的點，并且計算最佳的剛性變換，使得下述距離最?。?/p>

其中，分別為第k步迭代的點集p，q中的對應點，t＝(r；t)，r∈so(3)，t∈r³為剛性變換，γ代表所有剛性變換的集合；其中so(3)代表三階正交群，||·||代表歐氏度量，符號代表函數復合。(所有參數及符號定義下同)

具體來說icp主要包括下述步驟：

(1)對每一個pi^k，計算參考點集q中對應點使最?。?/p>

(2)計算旋轉矩陣r^k與平移向量t^k，極小化距離：

(3)更新第k+1步源點云中點的位置：

(4)計算l²距離誤差：

(5)如果d^k+1不小于給定的閾值δ，返回到(1)執(zhí)行，直到d^k+1＜δ或迭代次數大于預設的最大迭代次數為止；

最終，得到的變換為：其中ti為第i次迭代中求得的剛性變換。

因點集p為肋骨骨架內的骨骼點v^∞，而點集q為與之對應的肋骨二值分割結果(看作點集)，通過該算法，即實現(xiàn)將肋骨骨架與分割結果對齊。

利用cpr方法展開肋骨：

(1)應用得到的肋骨骨架(即通過應用最近點迭代矯正后的骨架)，均勻選取上面的3d采樣點；

(2)在每個采樣點上計算基于相鄰點和前一點的正交觀測方向，在每個點附近采樣一個平面；

(3)在每個點上重復上述過程，取一個切片，重復第二步；

(4)把采樣的切片相互堆疊在一起形成一個直卷；

(5)重復每根肋骨從而把肋骨展開。圖3和圖4是過程四得到的肋骨xz和yz方向的展開圖，可以從圖上觀察到骨斷裂處。

最后，需要注意的是，以上列舉的僅是本發(fā)明的具體實施例。顯然，本發(fā)明不限于以上實施例，還可以有很多變形。本領域的普通技術人員能從本發(fā)明公開的內容中直接導出或聯(lián)想到的所有變形，均應認為是本發(fā)明的保護范圍。