本發明涉及醫學圖像分割領域，尤其是在醫學圖像中對肺分割的方法及裝置。

背景技術：

肺分割是數字胸片圖像處理中的一個重要環節。分割結果的好壞直接影響后續對病灶的檢測及分析。胸片最大的缺陷在于組織重疊。人體結構中組織的密度越高，吸收X射線的能力就越強，因此胸部的肋骨吸收X射線就較多，在圖像中呈白影；由于肺部含氣體，所以密度較低，吸收X射線較少，因此圖像上呈黑影。正因如此，在分割肺邊界時，容易出現鋸齒狀的分割結果。而且通過觀察圖像可以發現，肺的下邊界處存在角點，這里也是分割的難點。

常用的肺分割方法大體可分為兩類：一類是基于規則級的肺分割，主要包括閾值分割法、區域生長法、邊緣檢測法、形態學濾波法等。由于X光圖像的成像效果的原因，這類方法不能對肺區域進行精細分割；另一類是基于像素分類的肺分割方法，主要有遺傳算法、神經網絡、模糊聚類等學習方法，對圖像中的像素進行分類。

還有一種是基于綜合知識的分割方法：主動形狀模型法(ASM)。這種方法首先需要人工標定肺的輪廓作為參考標準，通過訓練圖像樣本獲取特征點分布的統計信息，并且獲取特征點允許存在的變化方向，實現在目標圖像上尋找對應的特征點的位置。但是這種方法是在二維信息上進行查找邊界，很容易受到 肋骨邊界的影響，出現鋸齒狀的分割結果。

技術實現要素：

本發明所要解決的技術問題是提供一種能解決或降低前述問題的肺分割的方法及裝置。

本發明為解決上述技術問題而采用的技術方案是：一種醫學圖像中肺分割的方法，其特征在于包括以下步驟：

根據訓練集中的M張胸片的肺輪廓，得到平均肺模板，其中M為大于或等于2的整數；

對擬被分割的肺部圖像預處理，獲得預處理后的肺部圖像；

提取預處理后的肺部圖像的肺邊界的二值圖像，并根據廣義霍夫變換進行初始定位，獲得相對應的霍夫定位位置；

將平均肺模板與霍夫定位位置對齊，獲得對齊結果；

利用動態規劃算法進行分割，并將分割結果反變換回原始坐標系中，完成肺區域分割。

進一步的，在M張胸片中所對應的每張胸片的肺輪廓上標記相同數量的特征點，且不同胸片上的特征點的分布一致。

進一步的，所述在每張胸片的肺輪廓的拐點處、最高點處及最低點處標記特征點。

進一步的，所述在每張胸片的肺輪廓上除拐點處、最高點處及最低點處之外的其他位置等距離標記特征點。

進一步的，所述平均肺模板按照以下步驟獲得：

a.從訓練集中的M張胸片中進行第一次選取，獲取第一肺區形狀L₀；

b.從訓練集中的尚未被選取的胸片中進行第二次選取，獲取N張肺區形狀 L₁、L₂…L_N；

c.使所述N張肺區形狀L₁、L₂…L_N經處理后得到相應的N張肺區形狀L₁₁、L₂₁…L_N1，分別與第一肺區形狀L₀對齊，計算得到肺區平均形狀L_ave，所述處理包括使N張肺區形狀L₁、L₂…L_N分別相對于第一肺區形狀L₀作旋轉、縮放、平移變換；

d.若c步驟的肺區平均形狀L_ave相較于第一肺區形狀L₀的旋轉、縮放、平移小于設定閾值，則將所述肺區平均形狀L_ave作為平均肺模板；否則，進入e步驟；

e.將第一肺區形狀L₀進行處理、所述N張肺區形狀L₁₁、L₂₁…L_N1進行處理后，分別與當前肺區平均形狀L_ave-i對齊，計算得到本次肺區平均形狀L_ave1-i，所述對第一肺區形狀處理包括：使第一肺區的形狀L₀分別相對于當前肺區平均形狀L_ave-i作旋轉、縮放、平移變換，得到相應的肺區形狀L_0i，以及使所述N張肺區形狀L₁₁、L₂₁…L_N1分別相對于當前肺區平均形狀L_ave-i作旋轉、縮放、平移變換，得到N張肺區形狀L_1i、L_2i…L_Ni；

f.若e步驟計算所得本次肺區平均形狀L_ave1-i與當前肺區平均形狀L_ave-i的旋轉、縮放、平移小于設定閾值，則將e步驟計算所得本次肺區平均形狀L_ave1-i作為平均肺模板；否則，將當前肺區平均形狀L_ave-i用本次肺區平均形狀L_ave1-i替代，將第一肺區的形狀L₀用肺區形狀L_0i替代，N張第二肺區形狀L₁₁、L₂₁…L_N1用N張肺區形狀L_1i、L_2i…L_Ni替代，且i＝i+1，其中i大于或等于2的整數，“＝”為賦值號，并返回步驟e。

進一步的，所述肺部圖像預處理包括以下方式：獲取原始圖像，對原始圖像做高斯濾波處理后得到濾波圖像，對原始圖像與濾波圖像做差，然后施加濾波圖像的灰度平均值，得到去背景圖像，然后對去背景圖像做雙邊濾波操作， 即得到預處理后的圖像。

進一步的，基于預處理后的圖像，建立與圖像邊緣輪廓相對應的若干個邊緣形狀模板；獲取每個邊緣形狀模板的重心點的坐標以及邊緣形狀上的點相對于重心點的坐標位置，應用霍夫變換法對邊緣形狀的二值圖像進行初始定位。

進一步的，應用邊界增強算子提取邊緣形狀的梯度圖像，再保留最大的15％的像素點生成邊緣形狀的二值圖像。

進一步的，使用邊界二值圖像上每一個不為零的點去匹配模板中的任一點，在霍夫空間中相對應的重心點的位置處加1；結束二值圖像中所有等于1的點的霍夫變換后，求出霍夫空間中每一點的累加和，并找到最大值的位置，該點就是使用當前模板定位的最準確的模板的重心位置；對每個模板的重心點的最大值做進行歸一化，然后確定最大值對應的模板位置即為定位位置。

本發明為解決上述技術問題而采用的技術方案是：一種醫學圖像中肺分割的裝置，其特征在于包括：

存儲單元，存儲有訓練集的數張胸片圖像及擬被分割胸片圖像；

顯示單元，顯示對應的胸片圖像；

輸入單元，對胸片中的肺部標記特征點，獲得相應的肺輪廓圖像；

圖像處理單元，其用于：

對肺輪廓圖像進行處理后獲得平均肺模板；

對擬被分割的胸部圖像預處理，獲得預處理后的肺部圖像；

提取預處理后的肺部圖像的肺邊界的二值圖像，并根據廣義霍夫變換進行初始定位，獲得相對應的霍夫定位位置；

將平均肺模板與霍夫定位位置對齊，獲得對齊結果；

利用動態規劃算法進行分割，并將分割結果反變換回原始坐標系中，完 成肺區域分割。

本發明對比現有技術有如下的有益效果：本發明通過霍夫變換進行擬被分割的肺區進行定位，并通過一個平均肺區形狀與定位位置對齊，得到肺部的初始輪廓，最大程度的縮小了初始位置與肺部實際位置的偏差，從而提高了分割的準確性；另外，本發明是一種自動化的肺分割的方法，工作效率高。

【附圖說明】

圖1為本發明實施例的肺分割的方法示意圖；

圖2為本發明實施例的全自動肺分割的具體流程圖；

圖3(a)為本發明實施例中一獲取平均肺模板方法的流程圖；

圖3(b)為本發明實施例中另一獲取平均肺模板方法的流程圖；

圖4為本發明實施例中用于肺部上邊緣的霍夫模板示意圖；

圖5為本發明實施例中霍夫定位示意圖；

圖6為本發明實施例中肺部對齊結果示意圖；

圖7為本發明實施例中肺部分段分割結果示意圖；

圖8為本發明實施例中右肺分割結果示意圖。

【具體實施方式】

下面結合附圖和實施例對本發明作進一步的描述。

請參閱圖1-3，本發明實施例的一種醫學圖像(例如DR圖像或X射線圖像)中肺分割的方法，包括以下步驟：

根據訓練集中的M張胸片的肺輪廓，得到平均肺模板，其中M為大于或 等于2的整數；

獲取擬被分割的肺部圖像；

對擬被分割的肺部圖像預處理，獲得預處理后的肺部圖像；

提取預處理后的肺部圖像的肺邊界的二值圖像，并根據廣義霍夫變換進行初始定位，獲得相對應的霍夫定位位置；

將平均肺模板與霍夫定位位置對齊，獲得霍夫定位結果；

利用動態規劃算法進行分割，并將分割結果反變換回原始坐標系中，完成肺區域分割。

實施例中的肺分割，是通過醫學處理裝置進行全自動的肺分割。所述訓練集中包括若干張胸片，胸片被存儲在醫學處理裝置的存儲單元內；所述訓練集包含100張胸片，每張胸片對應一對肺輪廓圖像。在其他實施例中，訓練集中胸片的數量M可以為80張、120張、150張等。經發明人多次實驗，認為100張胸片即可較好的滿足實際需求；當然，胸片的數量越多，所得到的平均肺模板則越好，但是在圖像分割中所需要的時間也越長。

選取M張胸片，在M張胸片中所對應的每張胸片的一個肺部(邊緣)上標記相同數量的特征點，且不同胸片上的特征點的分布一致。實施例中，在每張胸片上標記42個特征點，其中在每張胸片的肺輪廓的拐點處、最高點處及最低點處都需要標記特征點。在每張胸片的肺輪廓上除拐點處、最高點處及最低點處之外的其他位置等距離標記特征點。當然，特征點的個數也可以為36個、50個等。標記后，就得到了M個肺的輪廓(形狀)。但是手工標定是會存在一定誤差，而且這些模型會出現在不同位置，并具有不同大小不同旋轉角度。為了消除這些“非形狀”因素，求取平均肺模板，我們需要對訓練集進行對齊。

具體實施方式中，利用循環計算(比較)的方式獲得所述平均肺模板，以下為相應的兩種可選步驟：

第一種獲取平均肺模板的方式(步驟)：

a.從訓練集中的M張胸片中進行第一次選取，獲取第一肺區形狀L₀；

b.從訓練集中的尚未被選取的胸片中進行第二次選取，獲取N張(第二)肺區形狀L₁、L₂…L_N；

c.使所述N張肺區形狀L₁、L₂…L_N經處理后得到相應的N張肺區形狀L₁₁、L₂₁…L_N1，分別與第一肺區形狀L₀對齊，計算得到肺區平均形狀L_ave，所述處理包括使N張肺區形狀L₁、L₂…L_N分別相對于第一肺區形狀L₀作旋轉、縮放、平移變換；

d.若c步驟的肺區平均形狀L_ave相較于第一肺區形狀L₀的旋轉、縮放、平移小于設定閾值，則將所述肺區平均形狀L_ave作為平均肺模板；否則，進入e步驟；

e.將第一肺區形狀L₀進行處理、所述N張肺區形狀L₁₁、L₂₁…L_N1進行處理后，分別與當前肺區平均形狀L_ave-i對齊，計算得到本次肺區平均形狀L_ave1-i，所述對第一肺區形狀處理包括：使第一肺區的形狀L₀分別相對于當前肺區平均形狀L_ave-i作旋轉、縮放、平移變換，得到相應的肺區形狀L_0i，以及使所述N張肺區形狀L₁₁、L₂₁…L_N1分別相對于當前肺區平均形狀L_ave-i作旋轉、縮放、平移變換，得到N張肺區形狀L_1i、L_2i…L_Ni；

f.若e步驟計算所得本次肺區平均形狀L_ave1-i與當前肺區平均形狀L_ave-i的旋轉、縮放、平移小于設定閾值，則將e步驟計算所得本次肺區平均形狀L_ave1-i作為平均肺模板；否則，將當前肺區平均形狀L_ave-i用本次肺區平均形狀L_ave1-i 替代，將第一肺區的形狀L₀用肺區形狀L_0i替代，N張第二肺區形狀L₁₁、L₂₁…L_N1用N張肺區形狀L_1i、L_2i…L_Ni替代，且i-i+1，其中i大于或等于2的整數，“＝”為賦值號，并返回步驟e。

進一步的，若e步驟計算所得本次肺區平均形狀L_ave1-i與當前肺區平均形狀L_ave-i的旋轉、縮放、平移大于設定閾值，但是當i大于或等于90時，停止迭代計算，并以本次肺區平均形狀L_ave1-i作為平均肺模板。

以上N為小于或等于M-1的自然數，i為大于或等于1的自然數。

所述閾值為使兩個肺區形狀(圖像)對齊時，其中一個肺區的形狀相對于另一個肺區的形狀所做的幾何變換的值：平移距離小于或等于0.01像素，并且旋轉角度小于或等于0.001*pi/180，并且縮放尺度小于或等于0.001。

第二種獲取平均肺模板的方式(步驟)：

a.從訓練集中的M張胸片中做第一次選取，選取一張胸片，獲取第一肺區形狀；

b.從訓練集中未被選定的胸片中做第二次選取，選取至少一張胸片，獲取對應的(第二)肺區形狀；

c.通過旋轉、縮放、平移使所述第二次選取所得的肺區形狀分別與第一肺區形狀對齊，計算得到當前肺區平均形狀；

d.若c步驟的當前肺區平均形狀相較于第一肺區形狀的旋轉、縮放、平移小于設定閾值，則將肺區當前平均形狀作為平均肺模板；否則，進入e步驟；

e.從訓練集中未被選定的胸片中做第N次選取，選取至少一張胸片，獲得對應的肺區形狀，通過旋轉、縮放、平移將所述第N次選取所得的肺區形狀分別與當前肺區平均形狀對齊，并計算出至第N次已選取所有肺區形狀的平均形 狀，其中，N為大于或等于3的整數；

f.若至第N次已選取所有肺區形狀的平均形狀相較于當前肺區平均形狀的旋轉、縮放、平移小于設定閾值，則將至第N次已選取所有肺區形狀的平均形狀作為平均肺模板；否則，將至第N次已選取所有肺區形狀的平均形狀作作為當前肺區平均形狀，今N＝N+1，其中，“＝”為賦值號；并返回步驟e。

進一步的，若訓練集中所有胸片被選定完畢，則以最后一次所選胸片所得肺區形狀與當前肺區平均形狀計算所得肺區形狀作為平均肺模板。

所述閾值為使兩個肺區形狀(圖像)對齊時，其中一個肺區的形狀相對于另一個肺區的形狀所做的幾何變換的值：平移距離小于或等于0.01像素，并且旋轉角度小于或等于0.001*pi/180，并且縮放尺度小于或等于0.001。

通過以上方法計算出平均肺模板后，將其保存在存儲單元中，以便于后續使用。

所述肺部圖像預處理包括以下步驟：獲取被拍攝部位肺部的原始圖像，對原始圖像做高斯濾波處理后得到濾波圖像，對原始圖像與濾波圖像做差，然后施加濾波圖像的灰度平均值，得到去背景圖像，然后對去背景圖像做雙邊濾波操作，即得到預處理后的肺部圖像。

對預處理后的肺部圖像分別提取其上、下、外邊界的二值圖像，采用基于廣義霍夫變換進行初始定位。所述上、下、外(左、右)側邊界共同組成肺部邊緣形狀(肺部輪廓)。根據邊緣形狀，可建立若干個與之較相似的邊緣形狀模板。具體實施方式中，可建立3至8個邊緣形狀模板，優選的數量為5個邊緣形狀模板。具體的，可分別將一個邊緣形狀劃分為上邊界、下邊界、外側邊界等幾部分，并將這幾部分分別進行處理。

以上邊界定位為例：首先標記5個上邊緣形狀(上邊界)模板(參圖4)， 記錄每一個模板的重心點坐標以及模板上每一點相對于重心點的坐標位置。應用霍夫變換算法對上邊界二值圖像進行初始定位。霍夫變換算法，是對圖像進行某種形式的坐標變換，使得圖像上給定形狀曲線上的所有點經過變換以后都集中到霍夫空間的某些位置上形成峰值點，這樣就把對原圖中給定形狀曲線的檢測問題轉化為尋找霍夫空間中峰值點的問題。

以任一模板為例，應用邊界增強算子提取上邊界的梯度圖像，再保留最大的15％的像素點生成上邊界的二值圖像。使用上邊界的二值圖像上每一個不為零的點去匹配模板中的任一點，在霍夫空間中相對應的重心點的位置處加1。在霍夫空間中記錄的是可能匹配準確的模板的重心位置，因此在霍夫空間中累積的值越大，那么該點是圖像中形狀的重心點的可能性就越大。在空間中找到最大值點，我們就確定該點為形狀的重心點。結束二值圖像中所有等于1的點的霍夫變換后，求出霍夫空間中每一點的累加和，并找到最大值的位置，該點就是使用當前模板定位的最準確的模板的重心位置。重心位置確定后，模板的相應位置也就確定，這就是找到的上邊緣位置。

使用5個模板分別對肺區域上邊界二值圖像進行霍夫變換，在5個霍夫空間中會得到5個最大值，由于模板大小不同，需要對最大值進行歸一化。再在這5個最大值中找到最大值，再確定該最大值對應的模板，我們就將這個模板確定為最佳模板。

然后對外側、下邊界也是與上邊界同樣的操作進行定位，定位結果如圖5所示；即可獲得與肺部圖像相對應的霍夫定位位置。

將平均肺模板與霍夫定位位置對齊，獲得對齊結果(對齊圖像，參圖6所示)。如圖7所示，將對齊后的結果(共42個點)中的三部分：外側、下側及內側分別標記點，1-16點(外側)、17-24點(下側)、25-42點(內側)，對 上側我們使用的是霍夫定位后結果(也可對外側、下側也可直接利用霍夫定位后的結果)，將這四部分分別進行線性插值，然后在每一點的法線方向提取m個點像素值，組成法線矩陣，在法線矩陣上進行動態規劃分割，再將得到的分割結果反變換回原始的坐標系中，就完成了肺區域的分割。最后將這四部分的分割結果進行聯接并對邊緣進行平滑操作，就完成了肺的分割(如圖8所示)。

本專利申請中動態規劃分割的過程：

動態規劃算法中的路徑選擇依賴于當前狀態，同時也依賴于之前已選擇的狀態，其通常用來搜索目標的最優邊界。

在相應的霍夫模板上取每點的法線方向的像素，并組成一個新的矩陣，稱之為“法線矩陣”(類似于極坐標變換后得到的矩陣)，在生成“法線矩陣”的同時也存儲每一點在原始坐標系中的位置坐標，方便進行反變換。

根據對齊得到的位置，將對應的邊界進行二階曲線擬合，這樣可以比較容易的計算出曲線上每一點的法線。

在每一點的法線方向，以霍夫模板為基準，向上取m1個像素長，向下取m2個像素長，采集每點的像素值，組成一個m*n的法線矩陣，n為霍夫模板的長度，m＝m1+m2。對這個法線矩陣進行動態規劃分割。具體包括以下幾個主要步驟：

(1)獲取局部代價

在動態規劃算法中，局部代價由內部代價和外部代價組成。內部代價用來衡量邊界點的平滑程度，邊界越光滑，內部代價越小。外部代價用來衡量圖像梯度變化的大小，梯度越大，外部代價越小。

假設“法線矩陣”的大小為m×n，則我們將內部代價定義為：

E_int(i，j)＝|j-k|/(j+k)，j＝1...n，k＝1...n，i＝1...m (1)

其中，j和k分別為“法線矩陣”中第i列和i-1列上邊界點的縱坐標，并對內部代價進行歸一化。內部代價的大小表明了邊界的平滑程度。

外部代價使用法線矩陣的梯度圖像的相反數來表示：

E_ext(i，j)＝-G(i，j) (2)

所以總的局部代價是由內部代價與外部代價的加權和表示的：

E(i，j)＝ω_int×E_int(i，j)+ω_ext×E_ext(i，j) (3)

其中ω_int和ω_ext分別表示內部代價與外部代價的權重。

(2)計算累計代價

計算累計代價是一個動態累加的過程，每一列的累計代價都是前一列的累計代價與當前列的局部代價的累加和。

因為第一列不存在內部代價，所以在累計代價的第一列只有外部代價而沒有內部代價。之后，每一列都是迭代的計算累計代價。每一列的累加代價都是前一列在一定范圍內累計代價與當前點的局部代價和的最小值。K表示第i-1列的搜索范圍，使每相鄰兩列選擇的邊界點不會有太大的跳動。在計算每一點的累計代價的同時還要記錄一下當前點取最小值時的k的取值，這樣便于最優路徑的需找。

(3)背向搜索最優路徑

計算完所有的累計代價后，找到累計代價最后一列中的最小值作為初始點，根據記錄的k值向前進行搜索，找到最優路徑。再變換回原始坐標系中，就得到了分割結果。

本發明還提供了一種醫學圖像中肺分割的裝置，其特征在于包括：

存儲單元，存儲有訓練集的數張胸片圖像及擬被分割胸片圖像；

顯示單元，顯示對應的胸片圖像；

輸入單元，對胸片中的肺部標記特征點，獲得相應的肺輪廓圖像；

圖像處理單元，其用于：

對肺輪廓圖像進行處理后獲得平均肺模板；

對擬被分割的胸部圖像預處理，獲得預處理后的肺部圖像；

提取預處理后的肺部圖像的肺邊界的二值圖像，并根據廣義霍夫變換進行初始定位，獲得相對應的霍夫定位位置；

將平均肺模板與霍夫定位位置對齊，獲得對齊結果；

利用動態規劃算法進行分割，并將分割結果反變換回原始坐標系中，完成肺區域分割。

利用上文提出的方法和/或裝置，分別使用DR設備提供的23個數據、JSRT數據庫的247個數據進行測試。測試結果表明，霍夫變換對肺區域的上、下、外邊界的自動定位都比較準確的，而且圖像質量較好，分割結果也很理想。

JSRT數據庫的圖像質量較差，在定位邊界及分割時都會有些偏差，但影響不大，結果可以接受。

將所有的分割結果分為4個等級：“很好”，“好”，“可接受”，“差”。共270例數據，其中大約95％的圖像分割好或很好。

雖然本發明已以較佳實施例揭示如上，然其并非用以限定本發明，任何本領域技術人員，在不脫離本發明的精神和范圍內，當可作些許的修改和完善，因此本發明的保護范圍當以權利要求書所界定的為準。