專利名稱:數(shù)字醫(yī)學(xué)圖像中解剖實體的灰度值模型和/或幾何模型的構(gòu)造方法
技術(shù)領(lǐng)域:
本發(fā)明涉及從數(shù)字醫(yī)學(xué)圖像訓(xùn)練集中構(gòu)造解剖實體的光度(也稱為灰度值模型)和/或幾何模型(也稱為形狀模型)的方法。
這些模型可以應(yīng)用于基于模型的分割處理中以分割新采集圖像中的解剖實體。解剖測量可以基于該分割的特征點。
背景技術(shù):
在放射實踐中,頻繁地使用幾何測量輔助對異常的診斷。為了進行這些測量,必需將關(guān)鍵用戶點放置在圖像中例如顯示裝置上顯示的圖像中的其相應(yīng)解剖標(biāo)記位置上。諸如兩點之間距離或者直線之間夾角的測量是基于這些關(guān)鍵用戶點的位置。另外,可以通過評估整體幾何以確定正?;虍惓?,這涉及對完整形狀進行分析。因此,需要自動和客觀地提取內(nèi)嵌在放射圖像內(nèi)的定量信息。
這種頻繁進行測量的示例為計算胸腔RX圖像中的心胸比率(CTR)(圖6)。這個比率定義為心臟在頂點水平處的橫向直徑與胸腔的內(nèi)徑(ID)的比例,即CTR=(MLD+MRD)/ID。心臟的橫向直徑由心臟左側(cè)的最大橫向直徑(MLD)和心臟右側(cè)的最大橫向直徑(MRD)組成。顯然,這種定義要求放射科醫(yī)生沿左和右胸廓(ribcage boundary)的內(nèi)圖廓尋找,以定位計算內(nèi)徑ID所需的點對(point pair)。該點對必需位于胸腔的最大內(nèi)徑處。同樣,必需找到左和右心臟陰影邊緣,以定位計算MLD和MRD之和所需的點。更具體地,這些點位于相對于脊柱中線的最遠端。尋找邊緣的過程需要放射科醫(yī)生進行解剖分割并在被分割的解剖上定位這些點(這個示例中總共有四個點)。該分割步驟,對于CTR計算的情形,相當(dāng)于描繪肺野(lung field)。
數(shù)字圖像中的許多其它測量遵從相似的方法,涉及分割解剖器官或?qū)嶓w,其中在該解剖器官或?qū)嶓w上分割幾何特征點和測量對象被確定并最終進行測量(圖5)。
參考心胸比率計算的示例,為了自動定位所需的點,需要在胸部放射圖像上自動地分割肺野的方法。
根據(jù)應(yīng)用,可以通過許多種方式解決分割問題。分割策略從早些年的計算機視覺的低水平策略進展到最近的基于模型的策略。
低水平方法依靠局部圖像算子,將具有不同光度特征的像素分離并將具有相似局部光度特征的像素分組。這兩種的示例為邊緣檢測和區(qū)域生長。盡管這些低水平方法性能差,但這些方法在絕大多數(shù)商業(yè)圖像分析工具中非常流行。主要原因為這些低水平方法易于理解和實施。然而對于復(fù)雜圖像數(shù)據(jù),例如醫(yī)學(xué)圖像中存在的以及如前所述的胸腔圖像的內(nèi)容為代表的復(fù)雜圖像數(shù)據(jù),這些方法的有效性是有限的。
更成功的方法包括有關(guān)待分割的形狀以及有關(guān)圖像中對象的光度或灰度級外觀的先驗知識。這些方法(稱為基于模型的方法)經(jīng)常是基于模板匹配。例如通過相關(guān)或使用廣義霍夫變換技術(shù)來匹配模板。不幸的是,對于醫(yī)學(xué)圖像的情形,模板匹配可能不成功。這是因為解剖對象呈現(xiàn)的形狀和灰度級外觀具有很大的可變性。
基于由Kass等(M.Kass,A.Witkin,和D.Terzopoulos,Snakesactive contour models,Int.J.Computer Vision,1(4)321-331,1988)提出的活性輪廓(contour)以及基于水平集合(J.A.Sethian,Level set methods and fast matching methods,Cambridge Univ.Press,Cambridge,U.K.1999)的方法能應(yīng)付更大的形狀可變性,但仍然不適用于許多醫(yī)學(xué)分割任務(wù),因為這些方法包括很少關(guān)于待分割對象的先驗知識。
手工參數(shù)模型克服了這個問題,但這些模型僅限于單個用途。
鑒于這些缺點,顯然,需要一種通用分割方案,可以使用示例進行訓(xùn)練從而采集有關(guān)待分割對象的形狀以及圖像中該對象的灰度級外觀的知識。由Cootes和Taylor(T.F.Cootes,C.J.Taylor,D.Cooper,J.Graham,Active Shape Models-their training and applications,Computer Vision and Image Understanding,61(1)38-59,1995)提出的主動形狀模型(ASM)滿足這種分割方案定義。由標(biāo)記點向量的主分量給出該形狀模型?;叶燃壨庥^模型描述了中心位于與對象輪廓垂直的每個標(biāo)記處的輪廓的歸一化一階導(dǎo)數(shù)的統(tǒng)計。通過最小化該一階導(dǎo)數(shù)輪廓與該輪廓分布之間的馬氏(Mahalanobis)距離,可以在新圖像中找到標(biāo)記的位置。這種算法開始于初始化評估,并執(zhí)行擬合程序,該擬合程序是標(biāo)記位移和形狀模型擬合的更迭。已經(jīng)設(shè)計出許多相似的方法,所有這些方法都采用了三步程序。首先,這些方法都使用形狀模型,保證產(chǎn)生看似可信的結(jié)果。其次,這些方法使用灰度級外觀模型,以將對象放置在邊界周圍或?qū)ο髢?nèi)部的灰度級圖形和從訓(xùn)練示例預(yù)計的灰度級圖形相似的位置上。最終,該算法通過最小化一些成本函數(shù)而擬合該模型。
由主動形狀模型提出的方法仍然面臨諸多限制。
第一限制為需要初始評估,稱之為初始定位問題。在一些情況下,需要在這個框架內(nèi)進行廣泛的搜尋以獲得合適的初始輪廓。
第二限制為形狀模型和灰度級外觀模型的交替使用。首先,使用灰度級外觀模型更新該分割評估。其次,將該形狀模型擬合到被更新的輪廓。不幸的是,如果灰度級外觀模型錯誤地定位標(biāo)記,則該形狀模型會受到誤導(dǎo)。
使用主動形狀模型分割方案的一些實驗呈現(xiàn)另一個問題。如果特定標(biāo)記的灰度級外觀模型被用于尋找該標(biāo)記的更佳位置,要求真實標(biāo)記位置位于該算法所探索的區(qū)域內(nèi)。在一些情況下,該算法會在另一個標(biāo)記的區(qū)域內(nèi)搜尋標(biāo)記。這意味著在錯誤區(qū)域使用了錯誤的灰度級外觀模型,導(dǎo)致不良分割。
本發(fā)明的目標(biāo)是提供構(gòu)造數(shù)字醫(yī)學(xué)圖像中一個或多個解剖實體的光度和/或幾何模型的改進的方法。
發(fā)明內(nèi)容
通過如權(quán)利要求1的方法實現(xiàn)上述方面。
在從屬權(quán)利要求中陳述了本發(fā)明的優(yōu)選實施例的具體特征。
該灰度值模型包括灰度值的概率分布以及灰度值在分布于解剖輪廓上一組標(biāo)記處的多個多尺度的灰度值導(dǎo)數(shù)。
該形狀模型包括沿解剖實體輪廓的連續(xù)標(biāo)記之間的每個連接向量的概率分布。
在下文中,術(shù)語灰度值模型、灰度值外觀模型、灰度級模型、灰度級外觀模型以及光度模型被當(dāng)作同義詞使用。
同樣,術(shù)語形狀模型和幾何模型被當(dāng)作同義詞使用。
術(shù)語特征圖像和特征也被當(dāng)作同義詞使用。
本發(fā)明方法的實施例通常被實施成計算機程序產(chǎn)品的形式,當(dāng)在計算機上運行時用于執(zhí)行本發(fā)明的方法步驟。
該計算機程序產(chǎn)品通常被儲存于計算機可讀取載體介質(zhì)例如DVD或CD-ROM中。備選地,該計算機程序產(chǎn)品采取電信號的形式,并可以通過電子通信與用戶通信。
在EP A 1349098中公開了由雙向鏈接外部圖像模型和內(nèi)部信息模型組成的,通過將測量對象和實體分組成計算機化的測量方案,而實現(xiàn)數(shù)字采集醫(yī)學(xué)圖像中測量自動化的方法。
在根據(jù)EP A 1349098的測量對話期間,從計算機提取測量方案并激活該測量方案。在被激活的測量方案的指導(dǎo)下,隨后對顯示的圖像進行測量。
在這種計算機化方法中,大量幾何對象被映射到數(shù)字圖像中,其它測量對象和最終測量實體(例如距離和角度)都是基于這些幾何對象。基礎(chǔ)的幾何對象典型地為關(guān)鍵用戶點,這些關(guān)鍵用戶點定義了幾何測量所基于的其它幾何對象。
然而上述專利申請沒有公開如何可以自動地實現(xiàn)該映射而無需用戶定位和操縱關(guān)鍵測量點。根據(jù)本發(fā)明的方法可以用于實現(xiàn)該映射。
提出基于離散點的對象表示,用于描述醫(yī)學(xué)圖像中的解剖輪廓。
公開了用于產(chǎn)生模型的兩種策略以及相關(guān)系統(tǒng)以及將這些模型應(yīng)用于分割真實圖像。這些策略,盡管分解成具有相同目標(biāo)的多個結(jié)構(gòu)單元,但構(gòu)造和采用了不同的光度和幾何模型知識。
首先,通過訓(xùn)練圖像而構(gòu)造新的灰度級外觀模型,編碼在標(biāo)記位置的光度知識。這個步驟利用了每個標(biāo)記周圍被采樣鄰域內(nèi)的強度關(guān)聯(lián)。
其次,構(gòu)造形狀模型,編碼標(biāo)記之間的幾何知識。這個步驟利用被分割對象的標(biāo)記之間的空間關(guān)聯(lián)。
在進一步的分割步驟中,光度和幾何知識被共同地用于分割新圖像上的一個或多個解剖輪廓。
結(jié)果的分割可以在對圖像執(zhí)行幾何測量的過程中被用于推導(dǎo)預(yù)期測量點的位置。
可以以直接或者間接的方式定義測量點。
在直接方式中,標(biāo)記的位置被用作獨立的測量點。
在間接方式中,通過分割過程得到的標(biāo)記通過曲線被互連,且該曲線上的特征點被定義成測量點。
灰度級外觀模型每個模型化系統(tǒng)中的第一步驟將標(biāo)記的幾何位置限制成位于線性路徑上,該線性路徑垂直于位于標(biāo)記中心上的現(xiàn)有輪廓或稀疏柵格。圖像中標(biāo)記點處的光度圖形被輪廓向量捕獲,其中在基于多個尺度下提取的圖像導(dǎo)數(shù)的特征圖像中采樣該輪廓向量。該輪廓采樣可以沿圖像內(nèi)的線性路徑,或者圍繞標(biāo)記點的圓形路徑。
特征圖像當(dāng)計算圖像內(nèi)所有圖像點的光度圖形時,獲得特征圖像,該特征圖像的大小等于原始圖像。相鄰運算通常被用于計算這些特征圖像。特定點的特征圖像內(nèi)特征的數(shù)值可被排列成特征向量。特征圖像的計算通常包括兩個步驟(a)使用一組卷積核進行線性過濾的步驟,和(b)涉及結(jié)合計算局部圖像統(tǒng)計的非線性點運算的多個子步驟的后期處理步驟。在步驟(b)中,可以單獨地執(zhí)行這些子步驟之一。
不同濾波器和不同類型的后期處理導(dǎo)致不同類型的特征圖像。例如,Laws構(gòu)造25個二維卷積核,作為重組理想特征輪廓例如Level、Edge、Spot、Wave和Ripple的向量對的外積。后期處理步驟包括非線性開窗口(windowing),其中將絕對過濾數(shù)值在更大窗口內(nèi)求平均以獲得圖像紋理能量測量。Unser使用的濾波器組構(gòu)造成對應(yīng)于公知變換的向量的外積,該公知變換是例如離散正弦(DST)、離散余弦(DCT)、離散偶正弦(DEST)、離散實偶傅立葉(DREFT)、離散實奇傅立葉(DROFT)變換,并用于后期處理通道方差的計算。最終,使用被調(diào)諧成不同空間頻率和不同取向的組合的對稱或反對稱Gabor濾波器構(gòu)造Gabor濾波器組,使得所有輸出具有相等的頻率帶寬。在這種情況下,該非線性點運算步驟可包括設(shè)定閾值。
使用基于多分辨率小波濾波器、Rank值過濾器或自適應(yīng)濾波器的特征圖像,可以測量其它圖像特性。
特征圖像基于局部無序圖像(LOI)。醫(yī)學(xué)圖像中解剖標(biāo)記周圍的圖像結(jié)構(gòu)呈現(xiàn)典型的灰度級圖形。在本發(fā)明中,在如后文中勾畫出的N個數(shù)學(xué)特征中捕獲這個結(jié)構(gòu)。使用特征圖像建造灰度級外觀模型的概念是稍早由B.Van Ginneken等,Active shape model segmentation withoptimal features,IEEE Trans.On Medical Imaging,21(8)924-933,2002提出的,該模型是基于由J.J.Koenderink和A.J.Vandoorn,Thestructure of locally orderless images,Int.J.of ComputerVision,31(2)159-168,1999提出的所謂局部無序圖像(LOI)。
可以利用泰勒展開式,使用K階多項式近似在感興趣點即標(biāo)記附近的位置x0的灰度值函數(shù)f。各項的系數(shù)由在x0的導(dǎo)數(shù)f(i)給出f(x)≈Σi=0K1i!f(n)(x0)(x-x0)i]]>這個函數(shù)近似中的所有信息都包括在導(dǎo)數(shù)f(i)內(nèi)。類似地,可由包括原始圖像的導(dǎo)數(shù)的圖像濾波器組(L,Lx,Ly,Lxx,Lyy,Lxy,…)描述該圖像??梢允褂迷S多模糊或內(nèi)尺度σ進一步計算每個導(dǎo)數(shù)圖像(derivative image)。通過計算后期處理步驟,最后獲得這些特征圖像,作為每個導(dǎo)數(shù)圖像的局部圖像統(tǒng)計。該局部圖像統(tǒng)計包括每個位置x0附近寬度為α的窗口內(nèi)的許多時刻(moment)的局部直方圖。在具體實施例中,計算第一和第二直方圖時刻、一直到二階的所有導(dǎo)數(shù)(L,Lx,Ly,Lxx,Lyy,Lxy)以及5個內(nèi)尺度(σ=0.5、1、2、4、8像素),形成總共N=2×6×5=60個特征圖像。計算直方圖的窗口尺度是根據(jù)內(nèi)尺度,即,α=2σ。
在特征向量中直接在行和列中使用這些導(dǎo)數(shù)(即,x和y)的缺點為平移和旋轉(zhuǎn)缺乏不變性(即歐幾里得變換)。除非這些圖像抗旋轉(zhuǎn)成標(biāo)準(zhǔn)取向,否則這些算子只能在具有相同取向的示例上進行訓(xùn)練。笛卡兒微分不變式(CDI)描述了與所選擇的笛卡兒坐標(biāo)系無關(guān)的圖像微分結(jié)構(gòu),因此可以省略抗旋轉(zhuǎn)步驟。由高斯微分算子構(gòu)造CDI在L.M.J.Florack等,Scale and the differential structure of images,Image and Vision Computing,Vol.10(6)376-388,1992中被描述。
下面為使用一直到二階導(dǎo)數(shù)的CDI(L,Lxx+Lyy,Lx2+Ly2,Lx2Lxx+2LxyLxLy+Ly2Lyy,Lxx2+2Lxy2+Lyy2)。同樣地,在每個像素位置,在多個內(nèi)尺度σ使用這些算子并計算范圍為α的窗口內(nèi)局部無序直方圖的頭兩個時刻,從而構(gòu)造特征圖像,其中α和內(nèi)尺度相關(guān)聯(lián),例如α=2σ。這些算子與圖像幾何的歐幾里得變換無關(guān),但仍然與所選擇的尺度參數(shù)σ有關(guān)。
輪廓和特征相似性測量灰度級外觀模型是基于平均輪廓和方差-協(xié)方差矩陣而建立,并利用輪廓這些點之間的強度關(guān)聯(lián),這將在下面解釋。輪廓處的強度由前面所勾畫出的強度特征向量表征。
為了將一個輪廓和另一個輪廓比較,需要知道輪廓中每個點是如何與任何其它點關(guān)聯(lián)的,并測量該輪廓的特定點的數(shù)值相對于任一其它點偏離平均值的程度。通過協(xié)方差測量捕捉這個統(tǒng)計關(guān)系,總是在兩個點(或這些點處的特征)之間計算該協(xié)方差。如果測量某一個點和它自身的特征值之間的協(xié)方差,則可以獲得方差。對于總的點數(shù)為2k+1時,若輪廓在輪廓的任意一側(cè)包括k個點,則存在(2k+1)k個關(guān)系。來自輪廓的點對之間協(xié)方差告知點對的第一個數(shù)值如何相對于第二對發(fā)生變化。當(dāng)協(xié)方差大且為正值時,表示當(dāng)?shù)谝稽c的特征值增大時,第二點成比例地增大。當(dāng)協(xié)方差大且為負值時,表示當(dāng)?shù)谝稽c的特征值增大時,第二點成比例地減小。當(dāng)協(xié)方差為零時,意味著這些特征值之間不存在系統(tǒng)耦合,即,這些特征值是相互獨立的。所有方差-協(xié)方差對可以排列成協(xié)方差矩陣S,該矩陣沿主對角線對稱。協(xié)方差矩陣因此捕捉了輪廓中特征值對之間的所有結(jié)構(gòu)關(guān)聯(lián)。實踐中,協(xié)方差矩陣是由學(xué)習(xí)采樣構(gòu)建成的,這種情形下為許多標(biāo)記處的輪廓匯集以及訓(xùn)練圖像集中的許多強度特征。
協(xié)方差涉及如何在高維空間,即,給定輪廓中2k+1個點的空間內(nèi)計算距離,從而計算當(dāng)前輪廓和模型輪廓之間的距離。當(dāng)輪廓只包括值為x的單個點時,通過將當(dāng)前輪廓點和模型輪廓點的數(shù)值相減并取絕對值,即d=|x-x|,就可以將這些輪廓和模型點數(shù)值x進行相似性比較。當(dāng)輪廓包括兩個點時,二元向量x和模型二元向量x之間的距離可以計算成平面內(nèi)的歐幾里得距離,即d=(x-x‾)′(x-x‾)=|x1-x‾1|2+|x2-x‾2|2.]]>在所有方向上衡量歐幾里得距離都是等同的。然而,當(dāng)變量x1和x2具有不等的標(biāo)準(zhǔn)差σ1及σ2時,等距離線的輪廓則為主軸平行于坐標(biāo)軸的橢圓??梢允褂脴?biāo)準(zhǔn)差的倒數(shù)衡量這些差異,得到加權(quán)的歐幾里得距離dw=(x-x‾)′W(x-x‾)=|x1-x‾1σ1|2+|x2-x‾2σ2|2.]]>矩陣W包括變量的方差的導(dǎo)數(shù),即W=1/σ121/σ22.]]>在線性或圓形輪廓中相互靠近的像素的強度水平,當(dāng)其屬于具有平滑變化強度水平的相同解剖區(qū)域時,將是強關(guān)聯(lián)的。當(dāng)像素對被認為屬于不同組織時,關(guān)聯(lián)將成反比。在任一情形中,成像噪聲和不重要的解剖變化會引入像素自身強度水平的非零方差。
因此,如果方差也是關(guān)聯(lián)的并具有不同方差,則必需插入?yún)f(xié)方差矩陣的逆矩陣,得到所謂的馬氏距離dm=(x-x‾)′S-1(x-x‾).]]>馬氏距離衡量多維數(shù)據(jù)點x與其平均值x的距離,使得在相同的多變量正常密度輪廓上的觀察將具有相同的距離dm。這些馬氏等距離輪廓在二維空間內(nèi)呈橢圓的形式,在三維空間中呈橢圓體。顯然,輪廓2k+1中點的數(shù)目通常遠大于3,這種情況下等距離軌跡在多維空間中呈超橢圓體的形式。
灰度級外觀模型是從訓(xùn)練圖像集中獲得的,包括平均輪廓中的每個標(biāo)記以及每個特征的協(xié)方差矩陣。所儲存模型的總數(shù)因此正比于標(biāo)記的數(shù)目和特征的數(shù)目。
形狀模型在形狀模型化步驟中,考慮了沿輪廓的標(biāo)記之間的位置關(guān)聯(lián)。由于對具有相對確定性形狀的解剖結(jié)構(gòu)的分割是針對有關(guān)該形狀必須優(yōu)選編碼該確定性成分的知識,以及相關(guān)的解剖變形,這樣進一步優(yōu)選地忽略了不相關(guān)的形狀細節(jié)。
描述了兩種實施例來考慮位置關(guān)聯(lián)。第一個實施例是基于涉及同時使用所有標(biāo)記位置的整體法;第二個實施例采用局部法,模型化連續(xù)標(biāo)記的空間部署。
PCA分析和擬合在第一實施例中,由手動放置于排列訓(xùn)練圖像集中的標(biāo)記構(gòu)成標(biāo)記坐標(biāo)位置的方差-協(xié)方差矩陣。該方差-協(xié)方差矩陣中的每個對角元代表這些標(biāo)記的坐標(biāo)值的方差。非對角元代表標(biāo)記坐標(biāo)值對的協(xié)方差,即,表示兩個標(biāo)記的坐標(biāo)是如何彼此相對變化的。該變化可以系統(tǒng)地相對其每個平均值沿相同的“正方向”而形成大的正協(xié)方差,或者可以系統(tǒng)地相對其每個平均值沿相同的“負方向”而形成大的負協(xié)方差,或者相對于其每個平均值隨機地沿相對的方向而產(chǎn)生零協(xié)方差?,F(xiàn)有技術(shù)主成分分析(PCA)中已知的本征向量/本征值分析將檢測這些標(biāo)記相對于平均標(biāo)記位置的任何全局系統(tǒng)變化模式。PCA旨在通過只保留具有最大變化的坐標(biāo)分量而減小數(shù)據(jù)集合的維數(shù),因此PCA對原始坐標(biāo)位置進行統(tǒng)計學(xué)地去關(guān)聯(lián)(decorrelate)。通過投影到這些高度變化的子空間,可以由保持重要形狀特征的維數(shù)更少的系統(tǒng)近似這些相關(guān)的統(tǒng)計量。足夠數(shù)量的這些主要模式被保留并用于形狀模型,忽略了由于噪聲變化引起的其它模式。
標(biāo)記連接向量模型第二實施例旨在模型化在局部水平布置連續(xù)標(biāo)記。然而,這種約束并不意味著不捕捉全局形狀。給出起始點坐標(biāo)以及這些連接向量,從該起始點開始,通過迭代地將連接向量加到當(dāng)前的點,就可以完整地重構(gòu)形狀輪廓。在像素水平的形狀描述中,由例如鏈?zhǔn)酱a利用這種性能。該連接向量為該曲線的局部切向量的離散近似。當(dāng)考慮到連續(xù)標(biāo)記之間的方向向量時,即一階幾何導(dǎo)數(shù),可以描述對象的形狀而與對象的平移無關(guān)。通過考慮例如由兩個連續(xù)連接向量夾角表達的連續(xù)標(biāo)記間的曲率,即二階幾何微分,則可以描述對象的形狀而與該形狀的平移和旋轉(zhuǎn)無關(guān)(所謂的剛性體變換)。這種情況下,當(dāng)起始點和一階幾何導(dǎo)數(shù)(即,起始連接向量)已知時,通過迭代地累加連續(xù)切向量的差向量,可以完整地重構(gòu)該曲線?;谙袼厮降逆?zhǔn)酱a的曲率通常對噪聲是高度相關(guān)的,這是因為其尺度處于極精細的水平。這種缺點在本發(fā)明公開中得到解決,這是因為標(biāo)記間隔足夠遠,使得可以實現(xiàn)與噪聲無關(guān)并且精確地計算一階或二階導(dǎo)數(shù)信息。由于切向量或曲率向量的大小與尺度有關(guān),將訓(xùn)練集的每個曲線的大小規(guī)一化成例如單位元素(unity),可進一步將曲線描繪成相當(dāng)于歐幾里得相似變換。排列訓(xùn)練圖像集中手動放置的標(biāo)記之間的平移向量(即,一階導(dǎo)數(shù)或切向量)或平移向量的向量差(即,二階導(dǎo)數(shù)或曲率)的平均和協(xié)方差統(tǒng)計構(gòu)成了本實施例的形狀模型。
分割系統(tǒng)上述模型可以用于分割系統(tǒng)。
分割系統(tǒng)也可以根據(jù)分別包括具有相關(guān)目標(biāo)的構(gòu)建單元的兩種路徑進行組織。
第一構(gòu)建單元旨在將示例模型標(biāo)記位置的容許位置限制在目標(biāo)圖像中的多個最可能位置。
對于使用特征圖像中線性輪廓代表標(biāo)記的情形,可以采用最佳點表決方法,該方法在穿過示例標(biāo)記的垂直線上對候選位置集進行分級。
對于圓形輪廓的情形,可以在每個示例標(biāo)記和柵格交點位置考慮矩形搜尋柵格(search grid),并根據(jù)穿過所有特征圖像的圓形輪廓到模型圓形輪廓的組合馬氏距離對柵格交叉位置進行分級。
第二和第三構(gòu)建單元共同組成優(yōu)化步驟。使用與每個標(biāo)記的保留容許位置相關(guān)聯(lián)的成本測量,可以構(gòu)造成本矩陣。該成本包括與灰度級外觀中一致性相關(guān)聯(lián)的項??梢允褂眉訖?quán)PCA輪廓擬合在全局水平校驗與形狀模型的一致性。備選地,可以通過結(jié)合與標(biāo)記候選位置相關(guān)聯(lián)的形狀成本,以驗證該形狀一致性。動態(tài)編程算法可以用于構(gòu)造穿過該成本矩陣的最小成本路徑,進而得到路徑具有最低總成本的最終分割。
最后,最終分割上的特征點可間接地定義關(guān)鍵測量點的幾何位置,且可以基于這些關(guān)鍵測量點產(chǎn)生所有相關(guān)測量對象和測量數(shù)量。
備選地,分割中使用的標(biāo)記點可用于直接定義測量點,且可提供每個測量點的搜尋柵格以定位測量點的允許位置。
通過下述描述和附圖,本發(fā)明另外的具體實施例和相關(guān)優(yōu)點將變得顯而易見。
圖1為說明第一模型化系統(tǒng)(模型化系統(tǒng)I)的流程圖。
圖2為說明第一分割系統(tǒng)(分割系統(tǒng)I)的流程圖。
圖3為說明第二模型化系統(tǒng)(模型化系統(tǒng)II)的流程圖。
圖4為說明第二分割系統(tǒng)(分割系統(tǒng)II)的流程圖。
圖5為說明利用分割方法(其中使用了根據(jù)本發(fā)明產(chǎn)生的模型),基于被分割的解剖對象上特征點的測量計算的流程圖。
圖6為用于根據(jù)胸部放射照片測量心胸比率的測量方案。
圖7(a)說明了胸部放射照片中左肺的手動描繪,在圖7(b)中使用標(biāo)記集{pi}i=1n表達該左肺,本示例中n=14。
圖8說明了使用n個標(biāo)記集p1=(x1,y1),…,pn=(xn,yn)對肺野的離散表達。該圖像尺寸被規(guī)一化在0…1之間。通過該離散標(biāo)記插入連續(xù)曲線。
圖9說明了分別在標(biāo)記附近的直線上和半徑為rc圓上的多個點的灰度值采樣的線性以及圓形灰度級輪廓。
圖10示出了包括14個標(biāo)記的左肺形狀上一些連接向量的實驗評估分布(a)v1,(b)v2,(c)v3。
圖11示出了胸部放射照片的第一時刻特征圖像,對應(yīng)于一直到二階的導(dǎo)數(shù)(L,Lx,Ly,Lxx,Lyy,Lxy)以及等于0.5、2和8像素的尺度σ。
圖12示出了胸部放射照片的第二時刻特征圖像,對應(yīng)于一直到二階的導(dǎo)數(shù)(L,Lx,Ly,Lxx,Lyy,Lxy)以及等于0.5、2和8像素的尺度σ。
圖13示出了將本發(fā)明應(yīng)用于胸部放射照片上肺野的分割的示例(a)手動分割左肺野和右肺野,(b)自動分割左肺野和右肺野。
圖14說明了心胸比率的計算(a)使用手動肺野分割得到CTRmax=0.47,(b)使用自動導(dǎo)出肺野分割得到CTRautomatic=0.45。
具體實施例方式
將參考具體應(yīng)用,即醫(yī)學(xué)圖像中肺野的分割,詳細地解釋本發(fā)明。
對象表示在下述本發(fā)明方法的具體實施例中,圖像中的解剖對象被數(shù)學(xué)地表示成位于封閉該對象輪廓上的固定數(shù)目的離散的標(biāo)號點,即p1=(x1,y1),…,pn=(xn,yn)。
輪廓{pi}i=1n從p1前進到pn并返回到p1。因此,可以通過離散形狀向量x=(x1,y1,…xn,yn)T捕獲該對象。選擇坐標(biāo)系統(tǒng),使得該圖像區(qū)域內(nèi)的所有點都落在域
×
內(nèi)(圖7)。
下述分割方案需要許多訓(xùn)練圖像,其中手動確定形狀向量x。一旦在該數(shù)據(jù)集上訓(xùn)練該算法,則該算法可以在新圖像內(nèi)產(chǎn)生該形狀向量。
可以由來自點{pi}i=1n的曲線插值重構(gòu)解剖對象的連續(xù)輪廓。一階插值使用線段連接連續(xù)的點,形成多邊形近似。高階插值可以用于獲得更平滑的輪廓。與階數(shù)無關(guān)地,該連續(xù)表示可以用于導(dǎo)出曲線的空間性能,例如沿曲線的切線方向,或者曲線上特定點處的法線方向。
模型化系統(tǒng)I(圖1)灰度級外觀模型灰度級外觀模型(也稱為光度模型)捕捉在所謂標(biāo)記(也稱為標(biāo)記點)附近范圍內(nèi)的空間強度分布。通過使用灰度級輪廓已經(jīng)實現(xiàn)了該灰度級外觀模型的空間分量。從其構(gòu)造可知,可以采用其它方案。
例如被用于精確地描述肺野輪廓的標(biāo)記的數(shù)目n可以僅為14,例如標(biāo)記1表示每個肺野的最高點,標(biāo)記7置于左心或右心陰影及其相應(yīng)半橫隔膜的截面處,標(biāo)記9表示肋膈角。其它標(biāo)記可以等距離地置于這些主要標(biāo)記之間。顯然,可以使用更多數(shù)目的標(biāo)記。
線性灰度級輪廓(圖9)在每個標(biāo)記,描述標(biāo)記周圍的典型圖像結(jié)構(gòu)的灰度級外觀模型垂直于輪廓地被采樣。在標(biāo)記的任何一側(cè),使用特定步進尺寸采樣k個像素,這給出了長度為2k+1的輪廓向量。這種采樣方案的優(yōu)點為能夠模型化在該標(biāo)記處的線性強度梯度。這些輪廓可以被歸一化,使得輪廓中元素的絕對值的總和為1。采樣的方向可以被計算為將標(biāo)記與其前一個標(biāo)記或后一個標(biāo)記連接而得到線段上的法線之間的平均法線方向,或者計算為連接該前一個標(biāo)記和后一個標(biāo)記的線段上的法線。
灰度級特征圖像標(biāo)記周圍呈現(xiàn)的圖像結(jié)構(gòu)被捕捉為N個數(shù)學(xué)特征,這里所解釋的數(shù)目更高。通過在每個位置x0周圍寬度為α的窗口內(nèi)計算導(dǎo)數(shù)圖像在多個時刻的局部直方圖,由此獲得這些特征圖像。在優(yōu)選實施例第一和第二直方圖時刻中,計算一直到二階的所有導(dǎo)數(shù)(L,Lx,Ly,Lxx,Lyy,Lxy)和5個內(nèi)尺度(σ=0.5、1、2、4、8像素),得到總數(shù)為N=2×6×5=60的特征圖像。計算直方圖的該窗口尺度是根據(jù)內(nèi)尺度,即α=2σ。特征圖像集可以被看作是灰度值函數(shù)在標(biāo)記附近的灰度值分解。
胸腔圖像內(nèi)的計算特征圖像的示例在圖11(原始圖像以及5個導(dǎo)數(shù)圖像在三個尺度的第一直方圖時刻)和圖12(原始圖像以及5個導(dǎo)數(shù)圖像在三個尺度的第二直方圖時刻)中示出。胸腔圖像通常在圖像中顯示成處于豎直位置,關(guān)于其它放射檢查的解剖對象可能呈現(xiàn)非標(biāo)準(zhǔn)的位置以及旋轉(zhuǎn),這種情況下可以采用笛卡兒微分不變式構(gòu)造特征圖像。
灰度級外觀模型訓(xùn)練集的所有圖像被用于建立每個特征和每個標(biāo)記的統(tǒng)計模型。將在特征j(例如2D圖像中其總數(shù)為60)中在標(biāo)記i(例如胸腔肺野的2D情形中其總數(shù)為14)采樣的歸一化輪廓示為向量gi,j,從訓(xùn)練圖像評估gi,j的概率分布,得到平均值gi,j和協(xié)方差Si,j。對于長度為2k+1的輪廓向量,該協(xié)方差矩陣的大小為(2k+1)×(2k+1)。特定標(biāo)記i的灰度級外觀模型包括平均輪廓gi,j,j=1…N和協(xié)方差矩陣Si,j,j=1…N??偡我暗幕叶燃壨庥^模型包括所有平均輪廓gi,j,i=1…n,j=1…N和協(xié)方差矩陣Si,j,i=1…n,j=1…N的集合。
灰度級成本在新特征圖像j中點p處為特定標(biāo)記i采樣的新輪廓gi,j之間的馬氏距離計算成hi,j(p)=(gi,j(p)-g‾i,j)TSi,j-1(gi,j(p)-g‾i,j)]]>馬氏距離越小意味著輪廓gi,j(p)源于平均值為gi,j且協(xié)方差為Si,j的高斯分布的可能性越大。因此,馬氏距離可用作灰度級成本測量,用hi,j(p)表示。該成本測量為給定圖像中位置p的函數(shù)。最可能根據(jù)特征j的標(biāo)記i的真實位置的位置為使hi,j(p)最小的點p。這樣就可以為每個特征j定義灰度級成本。
形狀模型由n個標(biāo)記點描述曲線輪廓,在肺野的情形中每個肺野具有一個曲線輪廓(圖8)。在s個訓(xùn)練圖像集合中手動地確定這些標(biāo)記點,得到一系列點(x1,y1),…,(xn,yn)。這些坐標(biāo)元組被依次置于表示該曲線的向量x=(x1,y1…,xn,yn)T中。接著,對這些訓(xùn)練圖像的形狀向量x進行主成分分析(PCA)。PCA將曲線投影到覆蓋了肺野的幾何變形的最重要模式的更低維空間內(nèi)??梢允褂?t<<2n近似每個曲線 x≈x+Φb其中x為平均形狀, 為與第t最大本征值相對應(yīng)的協(xié)方差矩陣x的本征向量。每個本征值決定該形狀模式的形狀的方差。由向量b表示的該曲線近似組成了該形狀模型,并將曲線x擬合到該形狀模型。這些本征形狀可以被看作是由標(biāo)記集表示的形狀的零階(位置的)幾何分解。形狀模型的作用是將對象的形變約束在由訓(xùn)練集設(shè)置的變化界限之間,例如相對于平均形狀的三個標(biāo)準(zhǔn)偏差。
分割系統(tǒng)I(圖2)分割算法最初將把代表肺野輪廓的曲線放置于圖像內(nèi),并通過迭代地更新曲線組成標(biāo)記的位置直至優(yōu)化標(biāo)準(zhǔn)達到最佳值為止而找到分割解決辦法。每個迭代包括下述步驟。
步驟1.最佳點表決假設(shè)具體標(biāo)記的當(dāng)前位置,將在垂直于輪廓的線上尋找該標(biāo)記的更佳定位。在該標(biāo)記的當(dāng)前位置的各側(cè)上評估了ns個位置。這些2ns+1個點之一將被表決成最佳下一個點。N個特征集合的每個特征對該2ns+1個點之一貢獻一個表決。該2ns+1點中被選擇的點根據(jù)馬氏距離標(biāo)準(zhǔn)具有最小的成本。對于每個該2ns+1點,由特征圖像j中被采樣的2k+1個點組成的輪廓被置于hi,j(p)方程中,且該方程中填充了恰當(dāng)?shù)钠骄鵪i,j和協(xié)方差Si,j。由該特征來表決具有最小馬氏距離的點。每個特征選擇一個點,導(dǎo)致在2ns+1個點即集合N1,…,N2ns+1(其中Σk=12ns+1Nk=N]]>)上劃分N個表決。
顯然,特征的數(shù)目N必需遠大于可選擇的點的數(shù)目2ns+1,否則可能分配到可選擇點集的點的表決太少。
特定標(biāo)記i所考慮的2ns+1個點的集合中的每個點接收到的表決數(shù)Ni,j,可以被看作是該具體點正是標(biāo)記所搜尋的點的置信等級。
步驟2.最小成本路徑優(yōu)化根據(jù)第一步驟,每個標(biāo)記可以分離地置換到2ns+1個可能的位置??蓪⒅眯诺燃塏i,j置于如下矩陣中R=N1,1...N1,2ns+1...Ni,j...Nn,1...Nn,2ns+1,]]>其中Σj=12ns+1Ni,j=N∀i=1,...,n.]]>更新輪廓涉及選擇矩陣C各行內(nèi)的一個元素;最小成本路徑(MCP)搜尋是指,將形成貫穿該矩陣的路徑的所有元素連接,使得成本標(biāo)準(zhǔn)最優(yōu)化(即最小化)。
一種合乎邏輯的選擇是取各行中置信最高的元素,這種途徑的缺點為它不考慮存在界外值(outlier)的情形。因此可能包括路徑內(nèi)的陡變。通過施加排除該陡變的平滑約束,這個缺點得到克服。
為此目的,首先如下所述地構(gòu)造成本矩陣C=1/N1,1m...1/N1,2ns+1m...1/Ni,jm...1/Nn,1m...1/Nn,2ns+1m,]]>其中,分母中的冪m為影響計算路徑的正數(shù)(m越大則對點的成本貢獻越小,因此指導(dǎo)路徑朝向該點)。其次,首先如下所述地構(gòu)造成本函數(shù)J(j1,j2,…,jn)=Σi=1nCi,ji+αΣi=1n|δi,ji-δi-1,ji-1|,]]>其中δi,ji為點i朝其目標(biāo)點ji沿該形狀法線的位移,α為加權(quán)因子。該方程中的第二項為平滑約束,該約束處罰如下情形的目標(biāo)點配置一個目標(biāo)點相對于(w.r.t.)前一目標(biāo)點的位移沿該輪廓陡變的位移。使用例如D.H.Ballard,C.M.Brown,Computer Vision,Englewood Cliffs,Prentice Hall Inc.1982,pp137-143中所描述的現(xiàn)有技術(shù)中已知的動態(tài)編程技術(shù),優(yōu)化成本函數(shù)J,得到使J最小化的目標(biāo)點(j1*,…,jn*)作為最低成本路徑(MCP)。
作為矩陣數(shù)據(jù)結(jié)構(gòu)的備選,還可以使用分層曲線圖。各層代表在標(biāo)記的垂直采樣圖像點陣列。使用一組弧從最后一層返回到第一層,由此增加該曲線圖。動態(tài)編程技術(shù)提取邊界作為穿過該曲線圖的最佳閉合路徑。
步驟3.未加權(quán)和加權(quán)的PCA輪廓擬合在上一步驟中,每個標(biāo)記i被賦予新位置ji*。序列(j1*,…,jn*)表示通過將其替換成x并通過求解Φb=x-x導(dǎo)出b而可以被擬合到形狀模型x≈x+Φb的輪廓。在該形狀模型中擬合x是指找到b值,使得x+Φb在高置信標(biāo)記中可以小誤差地近似x。
未加權(quán)的PCA輪廓擬合由于 因此方程數(shù)目多于未知量,該方程組為超定方程組,沒有精確解。因此,‖Φb-(x-x)‖p必需被最小化,以恰當(dāng)?shù)剡x擇p。不同的規(guī)范得到不同的最佳解。具體地p=2的情形,即最小二乘方問題,更容易處理,在現(xiàn)有技術(shù)中采用基于標(biāo)準(zhǔn)方程和QR因式分解的方法可以求解該方程。將實際點序列和擬合點序列之間的差異用擬合誤差e=x-(x+Φb)表示,則問題變?yōu)閷ふ易钚』痚Te的b。替換e,該優(yōu)化問題變?yōu)閙inbJ(b)=minbeTe=minb((x-x‾)T-bTΦT)((x-x‾)-bΦ).]]>如果J(b)=0,則J(b)最小化。J(b)的梯度J(b)可以寫成J(b)=-2ΦTW2(x-x)+2ΦTW2Φb.
因此,由優(yōu)化方程得到已知的標(biāo)準(zhǔn)方程J(b)=0ΦTΦb=ΦT(x-x),其中求解上述方程得到bb=(ΦTΦ)-1ΦT(x-x)該方程將圖像空間的數(shù)據(jù)向量x投影到PCA模型空間并返回近似系數(shù)b。可以通過評估形狀模型x=x+Φb而將近似結(jié)果b投影回到圖像空間。
加權(quán)PCA輪廓擬合當(dāng)將輪廓擬合到形狀模型時,相應(yīng)數(shù)目的表決 可以用作加權(quán)。因此,PCA模型將趨于賦予高置信點高的優(yōu)先級。置信等級低的點對該模型產(chǎn)生的輪廓的影響減小。優(yōu)化問題現(xiàn)在則包括最小化加權(quán)二次誤差(We)T(We)。使用和未加權(quán)情形中采用標(biāo)準(zhǔn)方程的相似方法,該問題的解b可表達成b=(ΦTW2Φ)-1ΦTW2(x-x)計算該方程以將x投影在PCA空間內(nèi),且使用形狀模型x=x+Φb可以往回投影。
迭代步驟1至3,直到滿足停止判據(jù),例如在連續(xù)迭代之間標(biāo)記位置基本上不再改變時。
未加權(quán)和加權(quán)PCA擬合都得到圖像內(nèi)標(biāo)記點集x的最終分割結(jié)果??稍谶@些點之間插入連續(xù)樣條曲線以解析地表示該分割。
模型化系統(tǒng)II(圖3)灰度級外觀模型灰度級外觀模型捕捉標(biāo)記鄰域內(nèi)的空間強度分布。在本實施例中,空間模型偏離在標(biāo)記鄰域內(nèi)采樣的圓形輪廓。用于精確地描述肺野輪廓的標(biāo)記的數(shù)目n可以僅為14,例如標(biāo)記1表示每個肺野的最高點,標(biāo)記7置于左心或右心陰影及其相應(yīng)半橫隔膜的截面處,標(biāo)記9表示肋膈角。其它標(biāo)記可以等距離地置于這些主要標(biāo)記之間。顯然,在根據(jù)本發(fā)明的計算方案中還可以使用更多數(shù)目的標(biāo)記。
圓形輪廓在標(biāo)記鄰域中選擇點的一個備選方式為,在以標(biāo)記為中心并具有至少一個半徑的圓上對點進行采樣。圖9中給出了圓形采樣的示例。如果標(biāo)記位于(x0,y0),則在半徑為rc的如下各點處采樣該圖像的灰度值函數(shù)(x0+rccos2πnc(i-1),y0+rcsin2πnc(i-1))i=1,...,nc.]]>采樣被置于長度為nc的輪廓向量中。nc的適當(dāng)選擇為12、8、6、4和3(對應(yīng)于30、45、60、90和120度角間隔)??梢酝瑫r考慮在一系列半徑處的多個圓形子采樣。表達成與圖像尺寸之比的無量綱量分?jǐn)?shù)的適當(dāng)?shù)陌霃絩c為0.05和0.025。
和線性輪廓相比,該方案的優(yōu)點為在鄰域周圍捕捉2D結(jié)構(gòu)。和線性輪廓相比,圓形輪廓可以更好地模型化諸如肋膈角落點或左心陰影與左橫隔膜交點的特殊解剖標(biāo)記,因為這些標(biāo)記具有肺野輪廓的不連續(xù)一階導(dǎo)數(shù)。在這些特殊解剖標(biāo)記附近區(qū)域中沿肺野輪廓的標(biāo)記,如果其長度太長,則還可能被線性輪廓模糊地表示,因為該輪廓可能跨過并列的肺野邊界。
圓形輪廓可以被歸一化,使得輪廓中元素的絕對值的總和為1。子采樣的角度可以從無子采樣(得到原始分辨率的該圖像的再次采樣)變化到對于圓形輪廓的情形,在主軸上例如僅僅4點的子采樣。
灰度級特征圖像標(biāo)記周圍的圖像結(jié)構(gòu)被捕捉為N個數(shù)學(xué)特征,正如本發(fā)明的發(fā)明內(nèi)容中所述。通過在每個位置x0周圍寬度為α的窗口內(nèi)計算導(dǎo)數(shù)圖像在多個時刻的局部直方圖,由此獲得這些特征圖像。在優(yōu)選實施例的第一和第二直方圖時刻中,計算一直到二階的所有導(dǎo)數(shù)(L,Lx,Ly,Lxx,Lyy,Lxy)和5個內(nèi)尺度(σ=0.5、1、2、4、8像素),得到總數(shù)為N=2×6×5=60的特征圖像。計算直方圖的該窗口尺度是根據(jù)內(nèi)尺度,即α=2σ。特征圖像集可以被看作是灰度值函數(shù)在標(biāo)記附近的灰度值分解。
灰度級成本在新特征圖像j中為特定標(biāo)記i采樣的新輪廓gi,j之間的馬氏距離計算成hi,j(p)=(gi,j(p)-g‾i,j)TSi,j-1(gi,j(p)-g‾i,j)]]>馬氏距離越小意味著輪廓gi,j(p)源于平均值為gi,j且協(xié)方差為Si,j的高斯分布的可能性越大。因此,馬氏距離可用作灰度級成本測量,用hi,j(p)表示。該成本測量為給定圖像中位置p的函數(shù)。根據(jù)特征j最可能為標(biāo)記i的真實位置的位置為使hi,j(p)最小的點p。這樣就可以為每個特征j定義該灰度級成本。
總體灰度級成本可通過組合所有特征j的所有灰度級成本的總和而得到,即hi(p)=Σj=1N(gi,j(p)-g‾i,j)TSi,j-1(gi,j(p)-g‾i,j)]]>該總和反應(yīng)了位置p處的灰度級圖形和標(biāo)記i處的預(yù)計灰度級圖形之間的相似性。根據(jù)灰度級外觀,最可能與標(biāo)記i一致的位置為下述優(yōu)化問題的解p0=argminρhi(p).]]>為了構(gòu)造灰度級外觀模型,必需從訓(xùn)練圖像中評估特征輪廓gi,j(對于所有標(biāo)記i=1,…,n,對于所有特征圖像j=1,…,N)的分布(gi,j,Si,j)。
形狀模型盡管灰度級成本驗證了灰度級圖形或特征的圖形,但仍定義了形狀成本以驗證兩個連續(xù)標(biāo)記之間的連接。
由n個標(biāo)記點描述曲線輪廓,每個肺野具有一個曲線輪廓。在s個訓(xùn)練圖像集合中手動地確定這些標(biāo)記點,得到一系列點(x1,y1)…(xn,yn)=(p1,…,pn)。這些坐標(biāo)元組被依次置于表示該形狀的向量x=(x1,y1…,xn,yn)T中。考慮連續(xù)標(biāo)記的評估位置對(pi,pi+1)。將形狀成本賦予向量v1=pi+1-pi,這反應(yīng)了vi的真實程度,即其概率分布。通過訓(xùn)練形狀評估v1,v2,…,vn的分布,假設(shè)在其平均值附近具有正態(tài)分布。vi的平均值和協(xié)方差分別用vi和Svi表示。在圖10中給出了這些向量分布的示例。這些向量分布可以被看作是由該標(biāo)記集代表的形狀的一階(正切)集合分解。
驗證兩個連續(xù)標(biāo)記pi,pi+1之間的連接向量vi(在平面內(nèi)可完全由其取向和長度表征的向量)的一種新穎形狀成本為vi與其平均值vi之間的馬氏距離fi(vi)=(vi-v‾i)TSvi-1(vi-v‾i).]]>由于具有大的馬氏距離而不可能發(fā)生的連接將具有高的形狀成本。另一方面,零成本將被賦予等于預(yù)期值的連接。
分割系統(tǒng)II(圖4)在訓(xùn)練階段期間采集有關(guān)圖像中對象的灰度級外觀的知識以及有關(guān)待分割對象的形狀的知識。這些知識隨后被用于根據(jù)下述步驟在新圖像中分割對象。
步驟1.矩形搜尋柵格構(gòu)造構(gòu)造每個標(biāo)記i的矩形搜尋柵格,以約束該標(biāo)記的可能位置。每個矩形柵格的特征在于幾何范圍以及柵格間隔的參數(shù)。
柵格間距rg應(yīng)該與灰度級外觀模型的圓形輪廓的半徑rc相對應(yīng)。適當(dāng)?shù)年P(guān)系為rg=fgrc,其中rc為半徑,fg為固定分?jǐn)?shù)。柵格間距的典型分?jǐn)?shù)為fg=0.5。
柵格范圍和柵格位置是針對每個標(biāo)記的,由xmin、xmax、ymin和ymax表示。真實(未知)標(biāo)記位置(x*,y*)被認為是位于下界和上界之間xmin<x*<xmaxymin<y*<ymax只有當(dāng)該柵格跨越整個圖像區(qū)域時才能滿足這些條件(即,xmin=y(tǒng)min=0以及xmax=y(tǒng)max=1)。更高效的方法是將搜尋約束在圖像的相關(guān)部分。假設(shè)該標(biāo)記的x坐標(biāo)的概率分布px(x)是平均值為x且標(biāo)準(zhǔn)偏差為σx(從訓(xùn)練形狀評估)的正態(tài)分布。在x方向上共同定義柵格范圍的柵格邊界xmin和xmax被選擇成使得∫xminxmaxpx(x)dx=fc]]>其中fc是代表條件xmin<x*<xmax為真的概率的分?jǐn)?shù)。分?jǐn)?shù)fc的恰當(dāng)值為0.995。如果施加的條件為該范圍圍繞平均值x是對稱的,即xmax-x=x-xmin,則上界xmax的值滿足1σx2π∫x‾xmaxexp(-(x-x‾)22σx2)dx=12fc.]]>下界則為xmin=2x-xmax。類似地獲得y坐標(biāo)的邊界。
步驟2.灰度級成本矩陣構(gòu)造灰度級外觀模型被用于尋找每個標(biāo)記的恰當(dāng)位置。如下所述地根據(jù)灰度級外觀模型選擇標(biāo)記i的m個最佳位置。
首先,根據(jù)前述步驟定義覆蓋圖像區(qū)域相關(guān)部分的矩形柵格,該圖像區(qū)域相關(guān)部分在每個標(biāo)記i的預(yù)計位置附近。
其次,對于每個柵格點,計算總體灰度級成本hi(p)。選擇和m個最低總體灰度級成本相對應(yīng)的點pi,1,pi,2,…,pi,m。
第三,構(gòu)造灰度級成本矩陣C,使得每行i包括標(biāo)記i的m個最可能位置的成本C=h1(p1,1)...h1(p1,k)...h1(p1,m).........hi(pi,1)...hi(pi,k)...hi(pi,m).........hn(pn,1)...hn(pn,k)...hn(pn,m).]]>每個標(biāo)記典型的最佳點數(shù)目m為10至40。
由于每個標(biāo)記的m個最佳點的選擇與鄰近標(biāo)記的m個最佳點集合相互獨立,因此會出現(xiàn)如下情況,即一個或多個這些點更靠近非鄰近的標(biāo)記。通過在下一步驟中考慮形狀成本,這些點可能在最后的輪廓分割中被忽略。
步驟3.最小成本路徑構(gòu)造確定分割對象的輪廓被簡化為通過選擇每行的一個元素在矩陣C中徹底地尋找路徑。將第i行中選定元素的指數(shù)寫成ki,則曲率變?yōu)辄c序列(p1,k1,p2,k2,…,pn,kn)。最佳路徑(k1*,k2*,…,kn*)是指使成本函數(shù)J(k1,…,kn)最小的路徑(k1*,k2*,...,kn*)=argmink1,...,knJ(k1,...,kn).]]>上面提出的模型允許進行大量的成本測量??紤]兩個連續(xù)標(biāo)記pi,ki和pi+1,ki+1,根據(jù)灰度級外觀模型的成本分量以及形狀模型的成本分量為-與標(biāo)記pi,ki及pi+1,ki+1相對應(yīng)的灰度級成本hi(pi,ki)和hi+1(pi+1,ki+1);-驗證從pi,ki到pi+1,ki+1連接的真實性的形狀成本fi(pi+1,ki+1-pi,ki)。
通過對總體灰度級成本和總體形狀成本的加權(quán)求和,將這兩種成本類型都結(jié)合到成本函數(shù)J(k1,…,kn)內(nèi)。加權(quán)因子為γ1的總體灰度級成本是標(biāo)記的單個灰度級成本之和hi(pi,ki),i=1,…,n。加權(quán)因子為γ2的總體形狀成本是各形狀成本之和fi(pi+1,ki+1-pi,ki),i=1,…,n。因此成本函數(shù)J(k1,…,kn)變?yōu)?
J(k1,...,kn)=γ1Σi=1nhi(pi,ki)+γ2Σi=1n-1fi(pi+1,k1+i-pi,k1)+γ2fn(p1,k1-pn,kn).]]>使J(k1,…,kn)最小的最佳路徑(k1*,k2*,…,kn*)稱為最小成本路徑。使用例如由G.Behiels等在Evaluation of image feature and searchstrategies for segmentation of bone structures using activeshape models,Medical Image Analysis,6(1)47-62,2002中提出的現(xiàn)有技術(shù)中已知的動態(tài)編程技術(shù),可以計算該最佳路徑。灰度級成本的典型加權(quán)因子為γ1=1,形狀成本的典型加權(quán)因子為γ2=0.25。
圖13示出了根據(jù)前述系統(tǒng)的三個胸腔圖像的自動分割,并將其與由經(jīng)驗的放射科醫(yī)生的手動分割進行比較。示出了每個肺野上的14個自動放置的標(biāo)記,同時還示出了穿過這些標(biāo)記的連續(xù)樣條曲線。
計算加速改進包括步驟1至3的分割方法可以被應(yīng)用于分割一個或多個解剖實體的輪廓的僅僅一部分。為此目的,只為這些連續(xù)標(biāo)記的子集構(gòu)造矩形搜尋柵格。該灰度級成本矩陣將包括數(shù)目與殘留標(biāo)記的數(shù)目相等的行。該最小成本路徑構(gòu)造將最小化成本函數(shù),該成本函數(shù)僅包括殘留標(biāo)記的灰度值項和形狀項。最終分割僅跟蹤被殘留標(biāo)記覆蓋的輪廓部分。當(dāng)一個標(biāo)記僅在部分解剖實體輪廓中感興趣時,或者當(dāng)需要定位僅部分特殊標(biāo)記以獲得具體測量時,部分分割的明顯優(yōu)點為計算速度。
加速分割計算的另一個算法改進是使用由粗到細的多分辨率方法?;诿總€標(biāo)記的多分辨率表示而構(gòu)造灰度級外觀模型和形狀模型,引入空間尺寸作為一個新的維數(shù)。為此目的采用了兩種備選,這兩種備選都在連續(xù)更精細分辨率水平上采用該分割概念。
在一個備選中,在分割期間,減少原始圖像中輪廓內(nèi)點的數(shù)目、搜尋柵格內(nèi)點的數(shù)目以及沿輪廓的標(biāo)記點的數(shù)目。通過在前一分辨率的附近范圍內(nèi)連續(xù)地提高搜尋柵格的分辨率(輪廓和搜尋柵格中所考慮的點的總數(shù)因此可以在每次迭代時保持不變)并采用更精細分辨率灰度值和形狀模型,于是可以提高這些標(biāo)記的位置精度。
在另一個備選中,可以使用多個水平的高斯和拉普拉斯金字塔及其相關(guān)導(dǎo)數(shù)特征圖像,實施該多分辨率方法。在粗糙水平上確定這些標(biāo)記的初始位置,并使用前一水平的中間位置作為下一水平的開始位置,由此跟蹤到最精細水平的最后位置。由于粗糙分辨率圖像包括更少的細節(jié),搜尋空間將可能包括更少的錯誤最小值,這使得可以在給定水平下實現(xiàn)更快的優(yōu)化并可朝其最后位置更快地定位這些標(biāo)記。
訓(xùn)練分割模型從圖像儲存器收集許多胸腔圖像,并將這些圖像依次表示在圖形用戶接口上。執(zhí)行下述步驟以建立分割模型—有經(jīng)驗的用戶通過描繪肺野的輪廓(使用鼠標(biāo)左擊)標(biāo)示出多個點而手動地分割每個被顯示圖像中的肺野。這些點無需沿肺野輪廓等間距地排布,唯一的要求就是這些點的集合整體上足夠高程度地近似該肺野。為了評估解剖擬合,將樣條曲線連續(xù)地插入穿過到目前為止已經(jīng)確定的手動定位點,直到通過鼠標(biāo)右擊使該曲線從最后一個點閉合到第一個點。通過將單個點朝新位置拖曳,可以進行手動分割調(diào)整。再次評估結(jié)果輪廓的解剖準(zhǔn)確性。
—接著,許多標(biāo)記將被自動放置到手動分割的肺野輪廓上。為了實現(xiàn)訓(xùn)練集的所有圖像上的相同標(biāo)記彼此映射,用戶定位幾個容易識別的標(biāo)記,通過等間距地放置于肺野輪廓上而獲得其它標(biāo)記。對于肺野輪廓的情形,多個容易識別的標(biāo)記為肺野最頂點、指示肋膈角的點以及心臟陰影和半橫隔膜交叉點,總共三個。接著選擇全部標(biāo)記,適當(dāng)?shù)倪x擇范圍為例如14至40。對于14個標(biāo)記的情形,點p1、p7和p9代表三個固定的標(biāo)記,5個點平均地分布在p1和p7之間,1個點位于p7和p9之間,另外5個點位于p9和p1之間。分別對左肺野和右肺野執(zhí)行這個步驟。
—接著,詢問有關(guān)訓(xùn)練階段的參數(shù)(a)用于訓(xùn)練(和分割)的圖像尺寸,典型值為256或512;(b)將由主成分分析(PCA)解釋的形狀變化的分?jǐn)?shù),典型值為0.9999;(c)線性或圓形輪廓中點的數(shù)目,圓形輪廓中的典型值為12、8、6、4或3;(d)作為圖像尺寸分?jǐn)?shù)的圓形輪廓半徑,典型值為0.04、0.03和0.02。
—訓(xùn)練灰度級外觀模型(a)在這些標(biāo)記附近分解該灰度值函數(shù),即計算N個特征圖像,例如如前所述的灰度值導(dǎo)數(shù)局部直方圖時刻,和(b)收集線性或圓形輪廓在標(biāo)記p1的灰度值分布,即計算gi,j,i=1…n,j=1…N以及協(xié)方差矩陣Si,j,i=1…n,j=1…N(即,對于特征圖像j內(nèi)的每個標(biāo)記i)。這個步驟產(chǎn)生統(tǒng)計灰度級外觀知識。
—訓(xùn)練形狀模型(a)計算包括在這些標(biāo)記中的零階(位置)信息的幾何分解,即計算平均形狀x以及在矩陣b中沿列方向排列的t個主要變形模式(本征形狀);(b)計算標(biāo)記連接序列中包括的向量分布(一階切線信息),即計算vi的平均值vi和協(xié)方差Svi。這個步驟產(chǎn)生統(tǒng)計形狀知識。
—儲存統(tǒng)計灰度級外觀和形狀知識,例如根據(jù)上面給出的基于模型的分割子系統(tǒng)來分割新圖像中的肺野。
將基于模型的2D圖像分割應(yīng)用于心胸比率(CTR)的計算自動肺野分割的一個用途為計算心胸比率(CTR)。CTR(圖6)定義為心臟的橫向直徑與胸腔的內(nèi)徑(ID)的比例CTR=MLD+MRDID]]>
MLD為心臟左側(cè)的最大橫向直徑,MRD為心臟右側(cè)的最大橫向直徑。這個比率是重要的臨床參數(shù),成人的變化范圍為39%至50%,平均值約為45%。高于50%的心胸比率被認為是異常的??赡艿脑驗樾牧λソ?、心包積液以及左或右心室肥厚。使用本發(fā)明中公開的自動肺野分割,可以自動地計算心胸比率。
參考圖8,通過選擇右肺分割的標(biāo)記p1和p9之間的擬合輪廓片斷上最靠左的笛卡兒點,從而獲得定義MRD的特征點;通過選擇左肺分割的標(biāo)記p1和p7之間的擬合輪廓片斷上最靠右的笛卡兒點,從而獲得定義MLD的特征點。將這些特征點的列坐標(biāo)相減并取絕對值,得到和MLD+MRD。類似地,通過下述步驟獲得ID選擇右肺分割的標(biāo)記p1和p7之間的擬合輪廓片斷上最靠左的笛卡兒點以及左肺分割的標(biāo)記p1和p9之間的擬合輪廓片斷上最靠右的笛卡兒點;將這些特征點的列坐標(biāo)相減;并取絕對值。
圖14示出了分割肺野的特征點的計算,以及心胸比率的計算(a)使用手動肺野分割得到CTRman=0.47,(b)使用根據(jù)本發(fā)明的自動導(dǎo)出肺野分割得到CTRautomatic=0.45。
將基于模型的分割和測量系統(tǒng)應(yīng)用于其它身體部分肺野分割中每個標(biāo)記的搜尋柵格的空間范圍是基于訓(xùn)練胸腔圖像集中所有相似標(biāo)記的位置而導(dǎo)出的。使用搜尋柵格用于約束給定標(biāo)記的候選位置的概念可以擴展到其它身體部分,這些身體部分在圖像中可能具有比胸腔圖像中肺野更大的位置、旋轉(zhuǎn)和尺寸變化,且和胸腔圖像中肺野相反地占據(jù)了整個圖像區(qū)域。
為了考慮到身體部分的這種位置、旋轉(zhuǎn)和尺寸變化,使用歐洲專利申請第04076454號標(biāo)題為“Method for automatically mapping ofgeometric objects in digital medical images”中公開方法的搜尋柵格定位點映射的概念可以和本發(fā)明結(jié)合使用。這種改進的方法于是對于分割整形外科放射照片中多骨的身體部分是特別讓人感興趣的,因為其特征在于在圖像中具有相對固定的形狀,這些部分可以被定位到圖像中非常明顯的多個標(biāo)記。
例如,在骨盆、髖部或整條腿檢查中,股骨輪廓可以被定位到熟知的標(biāo)記,例如稍大的轉(zhuǎn)節(jié)(trochanter)的尖端、膝蓋中心以及股骨頭部中心。這些定位點被用于建立模型定位點和實際圖像中選定的相關(guān)定位點之間的仿射變換。因為搜尋柵格包括收集排列在模型圖像中分割模型輪廓上標(biāo)記點附近矩形格子的候選點,通過將該變換應(yīng)用于模型點的坐標(biāo),每個組成的候選點依次被映射到實際圖像中。按照這個方式,在該點的最可能位置附近重構(gòu)特定標(biāo)記點的搜尋柵格。和前述方式相同地進行優(yōu)化過程,即通過優(yōu)化穿過候選位置選定組合的路徑的組合灰度值和形狀成本,每個搜尋柵格進行一次。具有最佳成本的路徑為該身體部分的最終分割。其它骨頭的示例為手和上肢骨頭、腳和下肢骨頭、骨盆和脊柱椎骨及結(jié)構(gòu)。
滿足這種類型的基于標(biāo)記分割的其它身體部分為MR圖像的3D CT中的軟組織器官。另一方面,這里可以采用基于逐個切片的方法,這種方法確定每個切片上解剖的分割輪廓,并將所有切片的結(jié)果組合成3D分割表面。另一方面,可以采用將本發(fā)明的概念擴展到3D體積的完全3D方法,該方法建立并應(yīng)用模型限制在3D搜尋柵格內(nèi)的3D標(biāo)記。腎、肺、心、肝、胃、脾、前列腺和大腦結(jié)構(gòu)為可以應(yīng)用本發(fā)明方法的器官示例。
在所有提及的應(yīng)用示例中,測量點可以是直接或非直接基于分割標(biāo)記的位置或者基于多個標(biāo)記的組合(例如伸長骨頭的皮層任意一側(cè)上兩個并列點的中點,一對該中點確定該骨頭的解剖軸),且可以將這些標(biāo)記點隨后輸入到諸如EP A 1349098中公開的測量系統(tǒng)內(nèi)。
權(quán)利要求
1.一種產(chǎn)生與數(shù)字醫(yī)學(xué)圖像中解剖實體相關(guān)聯(lián)的灰度值模型的方法,包括下述步驟主動在多個標(biāo)記點處采樣所述解剖實體的手動分割剖面;在每個所述標(biāo)記點的鄰域中采樣多個點;對于每個標(biāo)記點,將采樣點排列在輪廓標(biāo)記點周圍的鄰域內(nèi);計算至少一個特征圖像;計算每個標(biāo)記點和每個特征圖像的平均輪廓;計算每個標(biāo)記點和每個特征圖像的輪廓的協(xié)方差矩陣;確定用于所有特征圖像和所有標(biāo)記點的所述平均輪廓和協(xié)方差矩陣為所述解剖實體的灰度值模型,其中特征圖像為包括標(biāo)記位置周圍寬度為α的窗口內(nèi)所述圖像或?qū)?shù)圖像的局部直方圖的數(shù)學(xué)時刻的圖像表示。
2.一種產(chǎn)生數(shù)字醫(yī)學(xué)圖像中解剖實體的幾何模型的方法,包括下述步驟在多個標(biāo)記點處采樣所述解剖實體的手動分割剖面;計算連接向量或者連續(xù)標(biāo)記點之間的連接向量差;計算連續(xù)標(biāo)記點對的平均連接向量或平均連接向量差;計算連接向量或連接向量差的協(xié)方差矩陣;確定所述平均連接向量和協(xié)方差矩陣為所述解剖實體的幾何模型。
3.根據(jù)權(quán)利要求2的方法,其中所述灰度值模型和所述幾何模型是從圖像的多分辨率表示構(gòu)造而成的,并被應(yīng)用于所述圖像的多分辨率表示。
4.一種包含在計算機可讀介質(zhì)中的計算機程序產(chǎn)品,用于執(zhí)行權(quán)利要求1或權(quán)利要求2的方法的步驟。
5.一種計算機可讀介質(zhì),包括用于執(zhí)行權(quán)利要求1或權(quán)利要求2的步驟的計算機可執(zhí)行程序代碼。
全文摘要
產(chǎn)生一種灰度值模型,該灰度值模型編碼標(biāo)記位置的光度知識。這個步驟利用了標(biāo)記位置周圍被采樣的鄰域內(nèi)的強度關(guān)聯(lián)。產(chǎn)生一種幾何模型,該幾何模型編碼標(biāo)記之間的幾何知識。這個步驟利用了被分割解剖實體的標(biāo)記之間的空間關(guān)聯(lián)。
文檔編號A61B10/00GK1924929SQ200610126629
公開日2007年3月7日 申請日期2006年8月30日 優(yōu)先權(quán)日2005年8月30日
發(fā)明者P·德沃勒 申請人:愛克發(fā)-格法特公司