專利名稱：用于分割數字醫學圖像的方法
技術領域：
本發明涉及一種諸如乳房攝影照片之類的數字醫學圖像的上下文分割方法。
背景技術：
女性乳房的解剖結構隨著年齡而變化。在生殖年齡期間，乳房主要由導管、腺和脂肪組織組成。這散布有提供支撐和附著到胸壁的纖維組織。腺和纖維組織被共同稱為纖維腺組織。乳房腺組織被稱為乳房實質，并且由20..25個小葉(腺)組成，這些小葉負責生產乳汁，并且由集合在一起以形成較大導管的許多微管(導管)朝著乳頭引流。每個乳汁生產小葉包含一簇或一圈細胞。小葉和導管的部分由脂肪圍繞以進行保護，并且由纖維組織支撐。隨著年齡的增長，導管和腺單元經歷萎縮變化，并且逐漸由脂肪組織替代。乳房被將乳房組織附著到胸肌的韌帶固定就位。除了乳頭和圍繞它的乳暈之外，乳房到處被普通皮膚覆蓋。
纖維、導管和腺組織在X射線乳房攝影照片上呈現為暗的或“致密的”。另一方面，脂肪具有透明或明亮的外觀。術語“乳房攝影密度(MD)”和“乳房攝影模式”被廣泛地用來描述呈現在乳房攝影照片上的致密/明亮區域在乳房中的比例。在過去已經提出了乳房攝影實質模式的不同分級方法，例如Nottingham分級(5種模式，例如正常(N)，腺(G2，G1，G0)，發育異常(DS-輕度，DM-中度，DY-重度)，導管顯著(P1，P2)和不確定(IND))、Wolfe分級(4種類別)和Tabar-Dean分級(5種模式)。
乳癌的兩種主要類別是小葉和導管癌。
-小葉原位癌(LCIS)是在乳房的乳汁生產小葉中包含的細胞的數量、外觀和異常行為急劇增加的情況。術語“原位”是指癌癥的早期，并且被用來指示異常癌細胞存在但還未擴散超過它們初始發展處的組織的邊界。盡管并不認為LCIS是癌癥，但是診斷患有LCIS(也稱為小葉瘤形成)的婦女在一生中以后有發展成乳癌的更高風險。
-導管原位癌(DCIS)是乳房中早期癌發展的最常見情況。“原位”再次描述了還未移動到它初始發展處的身體的區域之外的癌癥。對于DCIS，癌細胞被限制在乳房中的乳汁導管，并且還未擴散到乳房脂肪組織中或身體的任何其他部分(例如淋巴結)。DCIS可以作為微鈣斑(稱為微鈣化)出現在乳房攝影照片上。
LCIS和DCIS都可以發展成侵襲性癌癥(浸潤性小葉癌或浸潤性導管癌)，其中癌癥擴散到乳房脂肪組織中或身體的任何其他部分(例如淋巴結)，這被稱為轉移。
到目前為止，乳房攝影法在乳癌的早期癥狀的檢測中已經成為最常用和最成功的工具，所述早期癥狀常常可以由微鈣化或腫塊的存在而得知。然而，由放射科醫生執行的視覺分析仍然是非常復雜的任務，并且開發了許多計算機輔助檢測/診斷(CAD)系統來支持他們的檢測和分級。實際上，CAD系統分別對有經驗的和無經驗的放射科醫生的檢測效率的影響例如在C.Balleyguier、K.Kinkel、J.Fermanian、S.Malan、G.Djen、P.Taourel、O.Helenon的“Computer-aided detection(cad)inmammographyDoes it help the junior or the senior radiologist？”，EuropeanJournal of Radiology 54(2005)(2005)90-96中進行了調查。在這兩種情況下，CAD系統已經證明是一種用于檢測的有效的支持工具，盡管其自主性仍有疑問。因此，由于問題的復雜性，自動或半自動系統仍然僅僅扮演放射科醫生的信令工具的角色。
在CAD環境中，圖像處理的作用之一將是檢測對于給定篩檢或診斷應用而言需要進一步處理的感興趣區域(ROI)。一旦檢測到ROI，后續任務將涉及所述區域的表征以及將它們分級成幾個類別中的一個。
乳房攝影病變的分級系統由美國放射學會提供，并且被稱為BI-RADS(乳房成像報告和數據系統)；除了X射線乳房攝影之外，包括關于乳房的超聲和磁共振成像(MRI)的部分。描述乳房攝影的發現的特征由描繪所述特征的線形圖來說明，之后是幾個乳房攝影的例子。
乳房攝影中ROI的例子是(a)鈣化、(b)腫瘤和腫塊、(c)胸肌、(d)乳房輪廓或皮膚-空氣邊界。分割是將圖像分成其組成部分、對象或ROI的過程。在圖像或其組成部分(即乳房攝影病變)的檢測、描述、識別或分級可以進行之前，分割是必需步驟。
諸如乳房攝影照片之類的放射圖像典型地由三個主要區域組成 -診斷區域包括對應于患者解剖結構(即乳房)的像素。通常，該成像區域的輪廓可以采用任何形狀。
-直接曝光區域是接收未衰減輻射的圖像區域。盡管該區域具有僅僅由噪聲破壞的恒定強度，但是入射能量(例如X射線源Heel效應)和接受器(例如計算機射線攝影中變化的存儲磷光體靈敏度)中的不均勻性會使該模式失真。在歐洲專利申請1256907中公開了一種方法，該方法從診斷圖像中回顧地估計這些整體不均勻性，并根據外推的背景信號來使所有圖像部分中的響應變平。
-如果有的話，準直區域作為高度衰減的像素出現在圖像上。這些區域的形狀典型地是直線的，但是圓形的或彎曲的準直形狀也可適用。
在放射圖像中的這些主要區域之間，可以考慮三個不同的區域過渡類型診斷/直接曝光、診斷/準直區域、以及直接曝光/準直區域的邊界。
分割算法旨在檢測并分離構成受分析的對象(一個或多個)的一組像素。這些技術可以根據應用于圖像的處理的類型被廣義地分類。基于區域的算法根據合適的相似性標準對圖像中的像素進行分組。在歐洲專利申請EP887769中公開了一種基于區域的算法，該算法通過根據灰度值直方圖的形心聚類對像素進行分組來分割直接曝光區域。基于邊緣的算法根據相鄰區域的灰度值差來分離圖像中的高對比度區域中的圖像像素。在歐洲專利申請610605和歐洲專利申請742536中公開了一種基于邊緣的算法，該算法檢測并描繪在單次或多次曝光圖像上的準直區域和診斷區域之間的邊界。在基于區域和基于邊緣的方法中，模型可以被用來限制所分割的圖像區域的外觀或形狀，以便遵守預定義的光度或幾何約束。該范例的例子是所謂的主動外觀模型和主動形狀模型(AAM和ASM)。
由于針對內乳房區域的分析，所以這些技術都不適用，原因是(a)它們是整體分割技術，其不適于乳房結構和病變的特定內容并且僅僅產生主要實體，例如乳房輪廓，或者(b)受到特定幾何或光度約束，其不適用于乳房結構和病變的乳房攝影外觀的大的可變性。所以，在當前的專利申請中，焦點集中于這樣的技術，所述技術可靠地分割乳房皮膚輪廓內部的結構和密度。
依賴于分割步驟的檢測步驟通常被認為是非常復雜的任務。例如，腫塊是比周圍組織更致密地成簇在一起的細胞群，并且可以在乳房攝影照片上由相對較小的強度變化來表示。此外，乳房攝影掃描記錄了乳房中存在的所有結構，其結構在尺寸、均勻性、位置和醫學意義上不同。最后，數字乳房攝影在其真正本性上固有的特征是由于現實世界和數字表示之間的轉換而引入的誤差，即量化噪聲。所有這三個因素增加了乳癌診斷任務的復雜性。
發明內容
本發明的一個目的是提供一種分割數字醫學圖像的方法，所述方法處理特別阻礙有用分割的數字乳房攝影照片的一般不均勻性。
上述目的由如權利要求1中所述的分割數字醫學圖像的方法來實現。在從屬權利要求中闡述了本發明的優選實施例的特定特征。
當使用像素灰度值作為局部圖像特征時本發明中的基本假設是，圖像由具有基本上一致的灰度值的區域組成。
在圖像不由基本上一致的灰度值的區域組成的情況下優選地使用多維特征，例如在局部導數的濾波器組的基礎上計算出的多維特征。
通常，聚類方案通過為每個像素一個分配表示級別的值來將輸入集(具有其強度值的圖像)分成一系列組，所述組被稱為簇。有效聚類的目的是加重存在于圖像中的對象之間的邊界，而不管尺寸、位置或可見度水平。聚類的最常用和最公知的技術之一是如普通教科書例如R.O.Duda，P.E Hart，D.G.Stork，Pattern classification，John Wiley and Sons，Inc.，2001中所述的K平均。該技術的特性和缺陷在乳房攝影分析的上下文中進行論述。
根據本發明，馬爾可夫隨機場模式被發展，并且作為擴展被應用于一般K平均方法。顯示出了簡單相關因子的引入如何顯著地改善結果。
在下面描述了K平均方案并且引入MRF算法。舉例說明了在乳房攝影數據分析上的主要應用。根據本發明的基于MRF的聚類模型的詳細定義與其參數和特征被一起公開。最后，由所述算法獲得的關于微鈣化和腫塊的檢測的結果被注釋，并且與經典K平均例程的結果進行比較。
根據后面的描述和附圖，本發明更多的優點和實施例將變得顯而易見。
圖1是來自MIAS數據庫的具有(a)微鈣化和(b)指示的腫塊的例子。
圖2是圖1(a)的放大片段，其中(a)發生微鈣化，(b)是其直方圖均衡形式。
圖3說明利用5個級別聚類的K平均。
圖4示出利用10個級別并且利用(左)和不利用(右)形心重算的聚類的結果。
圖5示出表示微鈣化的簇。
圖6是利用7個(左)和8個(右)級別聚類的結果；指示的區域是記錄的腫塊。
圖7是由低密度組織(左)和聚類的結果(右)表征的乳房的乳房攝影照片的片段。
圖8是由高密度組織(左)和聚類的結果(右)表征的乳房的乳房攝影照片的片段。
圖9示出n1鄰域系統和n1中的團集(clique)類型。
圖10示出n2鄰域系統和n2中的團集類型。
圖11是馬爾可夫隨機場聚類(MRFC)方法的操作的總框圖。
圖12是操作的詳細框圖。
圖13示出像素的3×3鄰域。
圖14說明MRFC方法的多分辨率逼近。
具體實施例方式 在當前的專利申請中，一種新的像素聚類模型被公開并且應用于數字乳房攝影照片的分析。所述聚類表示乳房攝影分析的更一般方法中的第一步，并且旨在創建簡潔數據集(簇)，以用于一方面是乳房攝影結構(例如胸肌和乳房輪廓)、另一方面是病變(例如腫塊和微鈣化)的自動檢測和分級。病變典型地是在乳癌的早期診斷中分析的第一癥狀。圖像像素由它們的強度(灰度級)或多維特征向量來描述。可以通過應用多尺度濾波器組或紋理算子來獲得多維分解。
引入了基于馬爾可夫隨機場(MRF)的技術，該技術適合用于執行由少量或有限的數據所表征的圖像的聚類。所提出的方法是統計分級模型，其基于圖像像素的統計和上下文信息的描述來標記圖像像素。除了評價源自K平均聚類方案的定義的像素統計之外，所述模型通過像素和它們的標記(上下文)之間的空間相關性的描述來擴展分析，因此相對于純K平均聚類的結果導致了分割輸出的不均勻性的減小。
乳房攝影數據 被稱為乳房攝影照片的乳房攝影檢查被用作一種篩檢工具，用于在無癥狀的婦女中檢測早期乳癌以及在有癥狀(例如硬塊、疼痛或乳頭溢液)的婦女中檢測和診斷乳房疾病。不能感覺但是可以在乳房攝影照片上看到并且告知可能的癌發展的異常之一是存在微鈣化，所述微鈣化是乳房組織中鈣的累積。微鈣化的尺寸從0.1mm到5mm之間變化，并且在乳房攝影照片上鈣群呈現為由較暗的不均勻背景圍繞的小點。告知可能的癌發展的異常的另一重要組是存在腫塊，所述腫塊是乳房組織中致密細胞的累積。在乳房攝影照片上細胞群呈現為由較暗的不均勻背景圍繞的亮區。
與所有X射線圖像一樣，乳房攝影圖像被各種性質的噪聲源破壞，例如(a)由于X射線光子通量的離散特性而造成的采集噪聲，(b)由于量化倉(bin)的有限范圍而造成的數字化噪聲，以及(c)由于投影圖像中解剖結構的疊加而造成的解剖學噪聲。因此，乳房攝影照片上的圖像區域本來是不均勻的，并且受到這些噪聲源的組合效應影響。
圖1(a)是來自公眾可得到的MIAS(乳房圖像分析協會)數據庫的在尺寸上減小的示例乳房攝影照片。所給出的圖像是發生微鈣化的情況；然而，異常的尺寸太小以至于不能容易地識別。在圖2(a)上更容易識別該鈣聯系(liaison)，在該圖中僅僅給出了相同乳房攝影照片的片段。仍然難以識別實際微鈣化和周圍組織之間的邊界。即使是相同片段的直方圖均衡形式的圖2(b)也不允許容易地分割所述結構。圖1(b)是發生腫決的情況，所述腫塊用所劃圈來指示。在識別實際腫塊和周圍組織之間的邊界上仍然有困難。
K平均統計數據聚類 數值數據的聚類形成許多分級和建模算法的基礎。聚類的目的是從大數據集中識別數據的自然分組以產生系統行為的簡潔表示。K平均例程根據關于實體所提供的數據將輸入集分成K簇。K平均聚類原則在諸如R.O.Duda，P.E.Hart，D.G.Stork的Pattern classification，John Wileyand Sons，Inc.，2001之類的關于模式識別的基礎教科書中進行了介紹。數據例如可以基于像素的灰度值通過特征空間中每個實體一個點來表示，或者基于多尺度濾波器組分解或一組紋理描述符，它可以是多維特征空間。開始，K級形心被定義，并且每個點被分配一個級別，其描述“最接近”所述形心。每次分配之后在屬于所述級別的點中重算形心。聚類以迭代方式重復所述分配，直到所有點和它們各自的級形心之間的距離之和被最小化。圖3給出了通過應用K平均方案并提供灰度級作為像素的描述來將圖像像素分成5個級別的結果。
可以清楚地看到，輸入圖像的空間不均勻性對所得到的簇集的形狀和不一致性具有強烈影響。K平均聚類例程的已知缺陷是，它不使用空間約束并且假設每個簇由恒定的強度來表征。后一種假設很少是正確的，并且也適用于乳房攝影圖像。減少問題的不完全解決方案可以通過將像素的坐標與灰度級值一起結合在所用數據中來實現。然而，在這樣的情況下，K平均將圖像分成K個空間相異簇，這是非常不希望的。也可以定義非常大數量的級別以滿足乳房攝影照片聚類的要求。然而，它將大大增加輸出的復雜性。
利用馬爾可夫隨機場建模 MRF理論為在貝葉斯框架內構造的視覺問題提供統計工具，其結合了結構和統計信息。MRF理論在S.Geman，D.Geman的公開文獻“Stochastic relaxation，Gibbs distribution，and the Bayesian restoration ofimages”，IEEE Trans.on Pattern Analysis and Machine Intelligence，Vol.6(1984)721-741中進行了介紹。
然而，為了與前述公布中的理論一致，聚類被定義為像素標記，像素被稱為位點(site)。
一組位點S＝{1，...，n}表示被標記的原始對象，即像素。向所有位點進行標記分配被定義為S上的隨機變量族的實現，并且被表示為{X＝x}＝{X1＝x1，...，Xn＝xn}(為了簡單起見，記作x＝{x1，...，xn})。描述位點的特征空間也被看作是具有其實現y＝{y1，...，yn}的隨機觀測場。通常，該組觀測向量不需要具有與該組位點相同的大小，例如當觀測被定義成用于成對位點時。然而，在本發明的上下文中，觀測的數量等于位點的數量。對于監督技術，位點的訓練集Str也被定義。在該情況下的訓練集由先前提取并標記(常常是人工的)的基元組成。本發明的技術是非監督技術。所以沒有定義用于聚類的訓練集。最后，標記集L={1，...m}被表示，其中每個標記對應于一個級別。值得一提的是，輸出可以由被分配相同標記的許多空間相異像素簇組成。
MRF中的最佳標記滿足最大后驗概率(MAP-MRF)準則，該準則要求標記的后驗概率P(X=x|Y=y)(在下文被表示為P(x|y))最大化，可以假設其遵循下列形式的吉布斯分布 其中配分函數Z是歸一化常數，U(x|y)是后驗能量，T是所謂的系統溫度并且在大多數視覺問題中被分配值1。吉布斯概率密度函數通過使用簡化的貝葉斯法則P(x|y)∝P(x)p(y|x)被結合在判定函數中，其導致下列形式 P(x)p(y|x)∝e-U(x)e-U(y|x)＝e-(U(x)+U(y|x))(2) 因此，來自等式1的后驗能量現在被表示為先驗能量項U(x)和似然能量項u(y|x)之和 U(x|y)=U(x)+U(y|x)(3) 吉布斯分布模型定義了位點i的標記，這取決于S中所有其他位點的標記。然而，被稱為Hammersley-Clifford定理的重要成果在吉布斯隨機場(GRF)和馬爾可夫隨機場(MRF)之間建立了一一對應等價。MRF模型是一種條件概率模型，其中位點標記的概率僅僅取決于給定鄰域內的位點標記。給定從本節開始時的記法，S的鄰域系統由N＝{Ni|i∈S)給出，其中Ni是與i相鄰但不包括f的所有位點的集合。
其中 依賴項的該“局部化”對于模型的復雜性的減小具有重要影響。相鄰的位點被分組成所謂的團集，取決于模型的復雜性，所述團集可以由任何數量的相鄰位點組成。團集的階指示所包括的位點的數量，并且它的增大通常增加了模型描述的復雜性以及其計算復雜性。
這里定義了一階的一元團集(其包括一個位點{i})和二階的二元團集(其包括一對相鄰位點{i，i′})。能量可以作為在團集上定義的局部勢能之和被估計。給定位點的標記，每個團集的團集勢能V(xi)和V(xi，xi′)和似然勢能V(yi|xi)被估計，其中i和i′標引S中的位點。
先驗能量是所有團集勢能之和 而在由標記xi調節的觀測是相互獨立的假設下，似然能量是所有似然勢能之和 通過最小化如等式(3)中所定義的能量來實現從等式(1)導出的MAP-MRF解 找到能量的全局最小值是一個組合問題，其通常使用以下算法來實現如J.Besag的“Spatial interaction and the statistical analysis of latticesystems”，J.Royal.Statistical Society Vol.36(1974)p.192-226中所介紹的確定性迭代條件模式算法(ICM)，或者如在S.Geman，D.Geman的“Stochastic relaxation，Gibbs distribution，and the Bayesian restoration ofimages”，IEEE Trans.on Pattern Analysis and Machine Intelligence Vol.6(1984)p.721-741中所介紹的隨機模擬退火(SA)算法。定義的組合特性暗示了進一步進行描述的模型的某些特征。
伊辛模型 伊辛模型是MRF中的自動邏輯(auto-logistic)模型的特殊情況，也被稱為多級邏輯(MLL)模型。它被表示為在相鄰位點上定義的簡單約束。當相鄰位點i和i′被分配兩個不同標記時的情況由團集勢能中的小常數來“處罰”，由此加強標記的均勻性；它可以以下列形式被表示 其中β是小正常數。伊辛模型的簡單性暗示了容易對標記過程包括某些約束。已知先驗信息可以在模型中被編碼以便增強所需的特性(一般(非)均勻性、定向(非)均勻性、紋理參數、...)。
操作的一般流程 圖11示出操作的一般流程的圖，其由下列步驟組成 -初始化和開始條件的評價。
初始化的過程由計算級形心的位置和它們的γ函數的相應值組成。
-在所有圖像像素上迭代 迭代的類型取決于在存在的幾種可能性中選擇的采樣方案。特別重要的兩種情況是 ○ICM-迭代條件模式。在該算法中沿著有序的像素列表來訪問像素。該列表從開始貫穿至結束并且從結束貫穿至開始。該方案中的標記也被排序，并且該算法依次評價它們中每一個的能量。在訪問所有像素之后，評價終止條件。
將ICM算法應用于乳房攝影分析的優點在于，一旦被正確地初始化，它將把較高標記分配給較亮區域，因而為分割區域的進一步分析自動地提供有用信息。
○SA-模擬退火。該算法在圖像中隨機選擇像素。然后還隨機選擇標記，并且僅僅為隨機選擇的標記和像素的當前標記評價能量。
模擬退火需要明顯更高數量的迭代。SA算法的優點歸因于這樣的事實，即它可以避免分割陷入局部最小值。缺陷可能是當應用于乳房攝影處理時SA會交換較高和較暗區域之間的標記，由此減少標記的含義。
-重算級參數。
每個標記改變觸發了兩個級形心和各自γ函數值的重算。要求重算的兩個級別的選擇由當前標記和新分配像素標記來確定。
-終止條件。
該算法評價在單次迭代期間改變其標記的像素的數量。重復迭代一直到不再有圖像像素改變其標記。乳房攝影照片的最終分割由對應于每個標記的圖像區域組成。
本發明的詳細描述 聚類方案的設計困難之一是級別的數量的選擇。一些例程被設計用于聚類成數量盡可能少的級別，并且具有最小的信息損失。然而在本情況中，盡可能少的方法將由于乳房攝影照片的特性而不會產生期望的結果。以非常簡單的方式，搜索乳房攝影照片上的腫塊的過程可以被描述為搜索很大網格的較亮區域。而且，被搜索對象的出現不明顯，并且通常它們的視覺特征非常類似于背景。實際上，強度上的微小差別是第一突出特征，并且眾所周知的是腫塊形成乳房組織內的稍微更亮的區域。在本情況下，數據由一組具有12位深度的強度圖像來表示，所以，數據的特征在于相對較低的分辨率。不過，已經設計了相干的和相對較簡單的MRF模型，其克服了所述問題并且執行有效聚類。
令S＝{1，...n}是位點即像素的集合，其中每個像素i由其灰度級值yi來描述。表示每個像素的標記分配的隨機變量的集合由x＝{x1，...，xn}表示，并且是所謂的配置。如前一節中所示，由于標記分配過程部分地取決于標記本身，因此配置x擴展了模型描述。最后，我們定義標記集L＝{1，...，m}，其中每個下標表示一個級別，該級別是一維灰度級特征空間中的子區域。
鄰域系統和團集 本發明的應用涉及在圖像像素的有限N1×N2矩形網格上定義為L＝{(m，n)∶1≤m≤N1，1≤n≤N2}的離散隨機場。當且僅當像素的鄰域ηmn僅僅包括作為(i，j)的鄰域的(k，l)的像素時，網格L的像素的子集的集合是L上的鄰域系統。這里使用鄰域系統的分級有序的序列。也被稱為一階鄰域系統的最近鄰域模型使用每個像素的最近四個鄰域，二階模型使用所有八個鄰域。與網格-鄰域對(L，η)相關聯的團集c被定義為L的子集，使得c由單個像素組成c＝{i}，或由一對相鄰位點組成c＝{i，i′}，或由三個相鄰位點組成c＝{i，i′，i″}組成，等等。這些鄰域系統和它們的相關團集在圖9和圖10中被描繪。
分割的特性可以與分配給每個團集類型的它們的控制參數一起被控制單團集位點的參數控制圖像中的每個區域類型中的像素的百分比；包含多于一個的位點的團集的參數控制聚類的大小和方向。
初始標記 基于MRF的例程的效率常常取決于起始點。在本申請中，該方案被定義為確定過程；所以，存在收斂到并不表示MAP-MRF解的局部最小值的危險。所以，通過為每個像素分配一個級別來執行初始標記，其灰度級最接近所述級別的形心。所得到的圖像是標記的散布集，在該圖像上容易看到當前噪聲和不均勻性的影響(參見圖3)。為了簡單起見，我們將一維特征空間分成等尺寸區域，并且將級形心置于所述區域的中心，即形心沿著特征空間均勻分散。初始標記根據預定數量的級別產生抖動圖像。在這一點上值得注意的是，初始標記等于K平均算法的第一次迭代的結果。
能量函數 在現有的問題中，能量函數扮演雙重角色。一方面，它必須引導例程給像素i分配標記，所述標記關于像素的強度值、因此其距級形心的距離而指向其級別。另一方面，它必須引導例程將相鄰像素分組成一致的、在空間上均勻的簇。
關于在2.3節中描述的和在等式(3)中以簡單形式給出的基本能量的形式，似然能量被定義為似然勢能之和。在本申請中，勢能表示像素距級形心的歐幾里德距離，從而產生 其中Cj，是級別j的形心的值。所定義的似然勢能導致標記的正確分配，不過，它由于在乳房攝影照片中存在的噪聲而導致空間非相干性。如果團集勢能未被定義(V(x)=0)，則所給出形式中的標記等效于K平均，并且歐幾里德距離被選作距離量度。由該例程所獲得的結果被顯示在圖3中。在等式(9)中定義的團集勢能表示能量函數的第二角色。該部分負責從標記中除去空間非相干性。它使用的唯一信息是一對相鄰像素的實際標記，并且采用如在S.Z.Li“Markov Random FieldModelling in Computer Vision”，Springer Verlag，1995中給出的改進的類伊辛形式 勢能“阻礙”了例程將不同標記分配給兩個相鄰像素。
能量函數和歸一化參數的一般形式 組合兩個所定義類型的勢能的問題起因于不同的值范圍。似然勢能的范圍取決于級別的數量、灰度級域的大小和級形心的位置，而團集勢能的范圍取決于鄰域系統的大小。所以，在試圖歸一化勢能的過程中，兩個參數α和β被引入以平衡(歸一化)能量函數。實際上，利用這兩個參數，可以在分配給勢能的類型的偏好方面調節模型的性能。通過如下對其擴展來修改來自等式(3)的能量形式 上面給出的形式假設參數α和β是常數。然而，在后面各節中顯示出，它們需要與級形心一起被重算。這兩個參數的函數在前面已被描述。它們對算法的影響是在期望性能的方面平衡例程并且導致它收斂。它們在根據兩種特殊情況的極端之間影響例程的性能 ·α＝0和β＞0 由于當使用各向同性鄰域(例如3×3鄰域)時β的值在所有方向上相等，因此標記的所得到的集合將僅僅由一個標記組成。在初始標記中出現最多的標記跨過網格分散。
·α＞0和β＝0 僅僅評價似然勢能，并且標記方案簡化為經典K平均算法。
值α和β對于能量函數的實際全局最小值具有直接影響。對于這兩個參數的人工設置可以迫使模型呈現某些特性(加強均勻性，最小化在像素值和它們各自的級形心之間的距離，...)。可以以這樣一種方式來設置參數，使得兩個相應勢能的“重要性”被平衡，即它們根據勢能的范圍被近似，以便在試圖歸一化勢能的過程中使二者相等。在一個級別內，即在分配給一個級別的像素的集合中，我們定義在灰度級域中的最大發生距離 如下計算所定義的最大可能距離α α(yi|xi)＝MaxDist(yi|xi)-1(14) 顯然，參數的該定義暗示了其對級形心的位置的依賴性，因此，當形心移動它們的位置時需要重算。以類似的方式計算參數β；僅僅團集勢能的域取決于假定的鄰域系統。由于我們將鄰域系統定義為8個鄰域，因此當具有標記xi的像素i被具有不同于xi的標記的8個像素圍繞時出現最大值。所以，我們計算 β＝8-1(15) 8個鄰域保證了在再分配標記時的各向同性；可以通過考慮如圖9和圖10中給出的定向團集來加入方向性。
改進的能量函數 如前一節中所示，參數β是常數，所以，它可以通過下列代入從等式中被除去 在等式(12)中，α也是常數，但是可能會顯示出，它的值遵循標記分配，而且，對于每個級別需要評價它。原因在于這樣的事實，即盡管級形心改變了它們的位置，但是級別的范圍也發生變化，這暗示了像素的強度和形心之間的距離分別遵循所述變化。由于α已被定義為級別內距離量度的歸一化參數，因此它必須與級形心一起被重算，并且由于新定義的參數γ是所分配標記和描述像素的特征的函數，因此能量函數采用改進的形式 級形心的值和γ′的值同時被重算；而且，遵循K平均的定義，在標記每次被分配給像素之后進行重算。
所述例程根據如在J.Besag的“Spatial interaction and the statisticalanalysis of lattice systems”，J.Royal Statistical Society，Vol.36(2974)p.192-226中所介紹的ICM采樣方案沿著圖像進行迭代。在每次迭代中，每個像素被訪問，并且對于S中可用的每個標記評價其能量。產生最低能量的標記如接下來所解釋的那樣進行選擇，各個級形心和函數γ(...)被重算，并且算法進入下一個像素。最后，每次迭代之后是計數改變其標記的像素，并且當計數結果低于預定義閾值時程序終止。
最低能量標記的選擇 通過省略等式(17)中的圖像環，由強度值y所表征的單個像素i的標記x的能量作為兩個部分之和被估計 其中V(yi|xi)是由函數γ(yi|xi)加權的似然勢能，并且V(xi，xi′)是團集勢能。Ni是如圖13中所繪的相鄰像素的下標的列表，代表性像素i用黑色繪出，并且其鄰域用灰色繪出。Ni是這些灰色像素的下標的列表。
似然勢能是在像素的強度值與由實際標記xi指示的形心的值之間的距離的量度。
團集勢能是所謂的罰值之和，并且當像素i和i′被分配不同標記時采用值一，否則采用零。這在當前的分割方法中是新穎的上下文術語，原因在于相鄰像素的標記影響當前像素的新標記的分配。它負責填充在分割標記的圖像中的局部間隙，由此減小由于噪聲或乳房組織的解剖學變化而引起的分段外觀。它因此克服了經典分割方法(例如K平均算法)的缺陷。
將所有可用標記依次分配給像素i并且估計相應的能量。選擇產生能量的最低值的標記，并且將該標記分配給像素。選擇“最佳”標記之后，該算法根據訪問圖像中的像素的順序進入下一個像素。
操作的詳細流程 圖12的詳圖擴展了圖11的總圖。
值“變化計數”形成算法的終止條件。在迭代期間當沒有圖像像素改變其標記時，算法終止穩定的標記。在每次迭代之前它被設置為零。
值“倒序”控制迭代的方向。在經典形式中，ICM采樣方案從第一個到最后一個、然后從最后一個到第一個地訪問有序像素的集合。一旦圖像像素已被“來回”訪問，迭代就被終止并且終止條件被檢查。
關于“變化計數”是否等于零的判定形成了算法的終止條件。
被稱為“選擇第一像素”的過程具有的含義是按順序選擇第一像素，即，當“倒序”被分配值零時選擇圖像的第一像素，而當“倒序”被設置為一時選擇圖像的最后一個像素。
過程“根據順序選擇下一個像素”也由“倒序”的值控制。在零的情況下，算法移到圖像的下一個像素，在一的情況下它進入前一個像素。
模型的性能 模型的性能在被應用于實際圖像上時關于兩個方面來分析(a)級形心重算的影響，以及(b)級參數的數量的選擇。在下文中使用具有微鈣化和腫塊的乳房攝影照片的代表性例子來評論這些方面。
形心重算。在每次標記分配之后，形心的重算在K平均方法的定義中具有其來源。如前所示，函數γ由于形心的重定位也需要重算。然而，當形心被設置成靜態并且重算未被執行時，最終結果(參見圖4)明顯發生變化。沒有重算的算法在較少數量的迭代中收斂，然而，結果指示了聚類的較低均勻性和所提取像素組的較低平滑度。
形心重算的原因是為了所公開算法的另一特征出現。在一些情況下，可能發生的是在標記的過程中級別中的一個將從圖像消失，當形心為靜態時這不會發生。實際上，由乳房攝影照片的高維度所產生的灰度級域的低密度和位點的高數量向我們保證，在初始標記之后每個級別將具有像素成員。然而，在進行中，如果算法檢測到級別還未獲得任何成員，即沒有像素被分配給所述級別，則它不能重算形心和函數γ(...)。在這樣的情況下，在初始標記之后兩個參數被重置為由計算產生的數值。重置這些值由這樣的事實證明，即如果一個級別從圖像中消失，則算法允許它再次重現。
級別的數量。被人工分配給算法的參數是級別的數量。實驗表明，小于10的數量對于聚類來說是滿意的，其將不會從最終結果中除去微鈣化和腫塊。然而，進行進一步的分析以評估描繪這些病變所需的級別的最小數量。
為此，使用10作為級別的數量，具有記錄的腫塊的圖像的集合被聚類。腫塊總是被正確地描繪，這允許計算它們的灰度級特性屬于腫塊的像素的灰度級范圍覆蓋由12位所表示的范圍的大約1/7.5。所以，確定了8是正確提取腫塊邊界所需的級別的最小數量。實際上，圖6示出分別使用7個和8個級別來聚類的結果。可以看到，聚類成7個級別并沒有提取可疑區域。然而，另一方面，8作為級別的數量卻成功地提取腫塊。
對于微鈣化適用類似的推理。同樣在這里，小于10個級別的數量對于聚類來說是滿意的，如果在圍繞圖像中期望微鈣化出現的位置的感興趣區域中局部地執行分析，其將不會從最終結果中除去微鈣化。圖5示出局部聚類成5個級別的結果。所指示的區域是保持它們的形狀和大小的所記錄的微鈣化。然而，當應用于整個圖像時，級別的數量必須增加以便保持這樣的小區域。而且，級別的數量的選擇由乳房攝影照片中呈現的乳房的密度來規定。由于X射線乳房攝影照片的疊加性質，因此不同組織彼此重疊，并且不同組織可以與微鈣化重疊。出現在乳房攝影照片中的微鈣化和其他結構之間的主要區別在于它們的大小和空間相干性微鈣化通常作為像素的圓形小組出現，其具有比周圍背景更高的強度值。
圖7和圖8給出了存在微鈣化的兩個乳房攝影照片的兩個片段，一個在主要由低密度組織組成的乳房中明顯存在鈣聯系，在另一個中高密度的乳房組織構成難以檢測的背景；微鈣化的兩種存在已由放射科醫生的檢查確認。當算法在整個圖像上運行之后提取片段。第二圖像明顯顯示出乳房攝影照片的更高密度；因此與用于第一乳房攝影照片的15個級別相對，它需要20個級別來提取微鈣化。
聚類的較高分辨率(即更多可能的級別)(a)導致在低可檢測條件中更好地檢測微鈣化，以及(b)導致更小微鈣化或對象的提取；因此它增加了算法的真陽性率和靈敏度。然而，級別的較低數量避免了提取不期望的對象；因此它導致低假陰性率并且增加了算法的特異性。在分析給出了由低組織密度表征的乳房的乳房攝影照片的情況下，在存在微鈣化時它產生更清楚的視圖；然而，如所示，它也在致密組織中導致較低的檢測率。
結果的分析讓我們定義了一個模型，該模型實際上是無參數的。實際上，8作為級別的數量在90％的情況下做得好。然而，這種類型的估計取決于數據集的特性，例如它取決于對于數字化膠片乳房攝影照片所用的掃描器的特性。盡管未根據經驗分配，但是級別的數量必須在掃描器數據或特定數據集的基礎上進行估計。
算法加速 可以通過使用圖像的多尺度分解來顯著地加速馬爾可夫隨機場聚類(MRFC)算法。在該分解中，高斯低通(LP)濾波被應用于原始圖像，并且該圖像隨后在每個維度上由因子2二次采樣(由符號↓表示)，從而產生第1級模糊圖像。該過程在每一等級i上被迭代許多次，從而產生所謂的高斯錐。在每次迭代i，上標度(由符號↑表示)粗糙細節的形式可以從前一精細細節的形式提取，從而產生帶通濾波圖像。這些帶通濾波圖像的序列被稱為拉普拉斯錐。通過將帶通濾波錐等級bi轉換成b′i并且重組它們，可以獲得增強的圖像。在圖14的虛線框中概括了這些操作。
從足夠粗糙的等級開始(例如圖14中的等級g4)，上面給出的馬爾可夫隨機場聚類算法被應用于尺寸減小的圖像上，從而產生分割的圖像s4。該等級的分割標記的初始化通過K平均聚類算法來獲得。由于在該等級上像素的數量被大大減小，因此獲得了ICM采樣方案的非常快的執行。
該等級的分割的圖像然后在每個維度上通過在每個維度上分割標記的像素復制由因子2擴展，并且該圖像構成下一個較高分辨率等級(其在圖1 4的例子中是等級3)的初始化。使用等級3的像素數據，MRFC算法確定新的分割標記，從而產生分割的圖像s3。由于已經從較粗糙等級正確地獲得了大多數標記，因此像素將僅僅局部地改變標記，并且因此該迭代也將在非常少量的時間內到達其終止條件。
在前一段中描述的該過程在下一個較精細等級上被重復，直到達到最高分辨率等級，從而產生原始圖像的分割s0。在圖14的虛線框的左邊描繪了這些操作的鏈。所述操作也可以應用于增強的圖像，其中對應于腫塊的微小區域關于周圍背景組織的對比度被增大。該操作鏈在圖14中的虛線框的左邊被描繪，并且在該例子中從被處理等級p4開始。通常，增強圖像的MRFC分割s′0將不正好等于原始圖像的分割s0。然而，在臨床環境中，原始粗數據大多數在PACS站不再可用，并且僅僅對比度增強的圖像將可用。對比度增強處理的例子是Agfa-Gevaert N.V.的MUSICATM算法。
其他分級數據結構可以被使用，例如四叉樹，正如在關于計算機視覺的普通教科書例如M.Sonka，V.Hlavac，R.Boyle的“Image Processing，Analysis，and Machine Vision”，2nd edition，Brooks/Cole PublishingCompany，1999中所述的。
可以通過使用快速多分辨率計算方案來緩解級別的合適數量的選擇，原因在于優化(例如ICM)是耗時的步驟。
為此，許多級別被輸入到加速算法，并且與每個級別數量相關聯的輸出分割被轉發到CAD算法以供進一步處理。
朝著在多維特征向量上操作的MRF的擴展 如先前概述的MRF模型由于根據所考慮的團集的相鄰像素的標記分配的相互作用而引入局部空間上下文。較大的空間上下文可以通過應用紋理算子來獲得，而不是僅僅使用像素的灰度值；因此，由于關于類似紋理模式而不是僅僅像素的強度的均勻性，所以可以推動朝著類似標記的分組。
在這樣的框架中，似然勢能現在是像素的多維紋理特征向量和實際標記xi所指示的多維級形心之間的距離的量度。團集勢能仍是所謂的罰值之和，當像素i和i′被分配不同的標記時它采用值一，否則采用零。
濾波器組可以由不同尺度和方向的一組導數濾波器組成；每個濾波器被稱為通道，并且提供多維特征向量的一個特征。例如，邊緣濾波器對應于一階導數；條形濾波器對應于二階導數濾波器。這兩個濾波器可以在多個尺度和方向上被計算。由于導數運算不是各向同性的，因此可以通過僅僅使用在所有方向上導致旋轉不變性的最大濾波器響應來使負責給定尺度的濾波器各向同性。由帶通濾波器組成的濾波器組的例子基于Gabor濾波器，該濾波器是高斯函數和復正弦的乘積。各向同性濾波器包括在不同尺度上計算的高斯濾波器和高斯拉普拉斯(LOG)濾波器。小波變換是可以作為數字濾波器組的樹實施的另一個例子。
當使用多維特征時，優選地使用在純歐幾里德距離上的其他距離量度來解釋給定級別的特征值中的可變性。
例如，可以使用加權歐幾里德距離，所述權重與標準偏差成反比例，或者可以使用馬哈拉諾比斯距離，其中權重由方差/協方差矩陣的逆矩陣表示。
對乳房攝影數據的應用 該算法的主要目的是描繪呈現在乳房攝影照片中的解剖學結構。該算法通過同時分析所有圖像像素的灰度級值來實現該目的，這允許在像素所屬的結構之間進行區分，同時增強被分配相同級別的標記的像素的空間分組。結果的特征在于解剖學結構的一致描繪，而不會使它們的邊界變形。
該算法的另一優點是當應用于原始輸入即未經處理的圖像時，實現了高性能聚類。實際上除了標記分組本身之外，以團集勢能的形式的依賴因子的公式提供了典型使用的預處理技術(即模糊化、直方圖銳化等)的替代，所述預處理技術被用來實現更一致或更平滑的輸入以供進一步處理。在乳房攝影照片分析中所述特征尤其重要，原因是乳房攝影照片的特征在于高灰度級分辨率的事實；所以任何類型的處理需要一種方法來處理被掃描的乳房攝影膠片或數字檢測器的自然噪聲。該算法減少了對預處理圖像的需要，并且以非常簡單的方式克服了圖像中存在的噪聲的問題。
該算法的又一個優點是一致數據集的提取以適合進一步的處理。該算法的輸出是在乳房攝影照片上的標記區域的集合。標記的值與圖像的灰度級范圍直接相關，即越高的標記是指越亮的區域。所以當涉及進一步處理時，比如腫塊檢測，根據標記的值來分析區域以提取ROI可能就已足夠，而不用大量分析原始圖像數據。
最后，該算法可以作為CAD(計算機輔助檢測/診斷)系統中的快速預處理步驟而被集成。模型的數學公式允許簡單的實施。當采樣方案被選擇成ICM時，該算法可以作為滑動窗口濾波器被實施。在ICM方案中，從開始直到末尾然后返回地訪問所有圖像像素。在每個點處，像素及其八個相鄰像素被用在計算中以考慮到相鄰標記的各向同性影響。所以，使用窗口寬度被設置為三的滑動窗口方案就已足夠，并且僅僅窗口函數需要被實施，所述函數將選擇中心像素的最低能量標記。
該方案提供了一種優化代碼以便實現所需的實時性能的容易的方式。
本發明特別適合于分割存在于乳房攝影照片中的乳房攝影組織類型，例如 -纖維腺盤 -脂肪組織 -皮膚 -胸肌 -病變的感興趣區域，例如包括微鈣化簇的區域或包括腫瘤塊的區域。
利用圖6作為例子，纖維腺盤由標記的小集合所表示的圖像區域的集合來表示。微鈣化或腫塊當存在時被表示為與小至僅僅一個分割標記相關聯的纖維腺盤中的疊加或嵌入區域。在乳房區域周圍的脂肪組織由與一個或兩個標記相關聯的細長區域來表示。胸肌在分割結果中自身呈現為對應于少量標記的圖像的角落中的區域，其中最小的密度標記具有典型為矩形的邊界。考慮到用于乳房攝影的CAD系統，這些區域的進一步自動分析現在在這些區域的幾何和放射測量特性的基礎上成為可能。
已經詳細描述了本發明的優選實施例，現在對于本領域技術人員而言顯而易見的將是，可以在其中進行大量修改而不脫離如由所附權利要求書中限定的本發明的范圍。
權利要求
1.一種用于分割數字醫學圖像的方法，包括以下步驟
a.在K平均聚類算法的基礎上計算每個像素的初始分割級標記，并且計算相關的級形心和伽馬函數；伽馬函數對能量函數中的似然勢能和團集勢能的重要性進行加權；
b.計算與依次將每個標記分配給所述像素相關聯的當前像素的能量量度，所述能量量度是以下兩個部分之和
-像素的特征值與由當前分割標記所指示的形心的值之間的距離的量度；
-被分配標記和至少一個相鄰像素的標記的差的量度；
c.為當前像素選擇產生最低能量的標記；
d.重算分割級形心和伽馬函數；
e.進入下一個像素，直到圖像中的所有像素都被訪問；
f.迭代步驟b-e，直到圖像像素的級標記不發生變化，以便產生圖像的最終分割。
2.根據權利要求1所述的方法，其中距離的量度是歐幾里德距離。
3.根據權利要求1所述的方法，其中所述距離的量度是加權距離。
4.根據權利要求1所述的方法，其中所述距離的量度是馬哈拉諾比斯距離。
5.根據權利要求1所述的方法，其中所述相鄰像素屬于各向同性鄰域。
6.根據權利要求5所述的方法，其中所述鄰域被限定為3×3鄰域。
7.根據權利要求1所述的方法，其中根據特定局部方向來選擇相鄰像素。
8.根據權利要求1所述的方法，其中像素由其灰度值表示。
9.根據權利要求1所述的方法，其中像素由多維紋理算子所提供的多維特征來表示。
10.根據權利要求9所述的方法，其中在局部導數的濾波器組的基礎上計算所述多維特征。
11.根據權利要求1所述的方法，其中步驟(d)和(e)的順序被交換。
12.一種用于分割數字醫學圖像的方法，包括以下步驟
a)計算數字圖像的多分辨率(MR)表示，
b)在所述MR表示中在粗分辨率等級選擇一個圖像作為原始圖像，并且將K平均聚類算法應用于所述原始圖像，
c)計算與依次將每個標記分配給所述像素相關聯的當前像素的能量量度，所述能量量度是以下兩個部分之和
-像素的特征值與由當前分割標記所指示的形心的值之間的距離的量度；
-被分配標記和至少一個相鄰像素的標記的差的量度；
d)為當前像素選擇產生最低能量的標記；
e)重算分割級形心和伽馬函數；
f)進入下一個像素，直到圖像中的所有像素都被訪問；
g)迭代步驟b-f，直到圖像像素的級標記不發生變化，以便產生圖像的最終分割；
j)將標記擴展到多分辨率表示的下一個更精細的分辨率等級；
h)迭代步驟(c)-(j)，直到原始圖像被分割。
13.根據權利要求12所述的方法，其中所述多分辨率表示是高斯錐表示。
14.根據權利要求12所述的方法，其中所述多分辨率表示是四叉樹。
15.根據權利要求12所述的方法，其中通過在所述多分辨率表示中的更高分辨率等級將標記局部復制到圖像來執行擴展。
16.根據權利要求12所述的方法，其中步驟(e)和(f)的順序被交換。
17.根據權利要求1所述的方法，其被應用于乳房攝影圖像。
18.根據權利要求1 2所述的方法，其被應用于乳房攝影圖像。
19.一種計算機程序產品，適于在計算機上運行時執行權利要求1所述的方法。
20.一種計算機程序產品，適于在計算機上運行時執行權利要求12所述的方法。
21.一種計算機可讀介質，包括適于執行權利要求1所述的步驟的計算機可執行程序代碼。
22.一種計算機可讀介質，包括適于執行權利要求12所述的步驟的計算機可執行程序代碼。
全文摘要
本發明公開了一種用于分割數字醫學圖像的方法。一種基于馬爾可夫隨機場(MRF)的技術被描述以用于執行由少量或有限數據表征的圖像的聚類。所提出的方法是統計分級模型，該模型基于圖像像素的統計和上下文信息的描述來標記圖像像素。除了評價源自K平均聚類方案的定義的像素統計之外，所述模型通過像素和它們的標記(上下文)之間的空間相關性的描述來擴展分析，因此相對于純K平均聚類的結果而導致分割輸出的不均勻性減小。
文檔編號G06T5/00GK101169868SQ20071018167
公開日2008年4月30日 申請日期2007年10月24日 優先權日2006年10月25日
發明者P·德維爾, G·貝希爾斯, E·尼森, R·德勒克, M·蘇利加 申請人:愛克發醫療保健公司