專利名稱：一種醫(yī)學圖像中腹部器官分割方法
技術領域：
本發(fā)明屬于醫(yī)學圖像處理及應用領域，涉及一種醫(yī)學圖像中腹部器官分割方法。
背景技術：
腹部器官的分割具有重要研究意義和臨床價值。首先，從各種類型的醫(yī)學圖像中將器官分割出來是進行可視化的第一步。比可視化更重要的，確定感興趣器官位置和區(qū)域在放療手術計劃中具有重要意義。只有掌握了器官的空間幾何形狀，才能制定出準確的放療方案和放療劑量。
Linda G. Shapiro等人1995年在《Pattern Recognition》雜志上提出利用器官形狀的先驗知識來動態(tài)調整灰度閾值進行分割，然后按照先易后難的順序逐個識別感興趣的器官。這種方法對每個器官采用的仍然是單一閾值分割，得到的區(qū)域往往是不連續(xù)的。John E.Koss等人1999年在《IEEE Transactions On Medical Imaging》雜志上提出選取圖像的紋理特征，采用Hopfield網絡進行聚類分割。這里分割得到的器官取決于Hopfield網絡中指定的聚類數(shù)目。Chien-Cheng Lee等人2003年在《IEEE Transactions On Information Technology In Biomedicine》雜志提出采用上下文相關神經網絡基礎上的模糊識別規(guī)則對各區(qū)域進行識別和合并。這樣一套模糊識別規(guī)則的建立，要求相當?shù)慕涷灪投啻卧囼烌炞C。Hyunjin Park等人2003年在《IEEETransactions On Medical Imaging》雜志提出則將關于腹部器官的解剖學形態(tài)圖譜通過配準集成到貝葉斯分割框架中[5]。該方法需要手工設置多個控制點，影響了方法的自動化。

發(fā)明內容
本發(fā)明的目的是為臨床應用提供一種在PC機上運行的腹部CT和核磁圖像的器官分割算法，其特點是同時進行圖像中腹部多器官自動化分割，實施簡單，無需人工干預。
算法的主要步驟包括步驟(1)，在PC機中，采用基于歸一化互信息的配準方法，將圖譜向通過CT或核磁掃描得到的個體圖像作整體配準，以消除圖譜和目標圖像間的整體差異；所述圖譜是指對已有的醫(yī)學圖象經解剖學專家手工分割后對每個器官區(qū)域做好標記的圖象，它可以為PC機提供人體解剖結構的參考信息；所述的歸一化互信息nMI(S，T*A)用下式表示
nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>其中，S為目標圖像，T*A代表對圖譜A做空間變換T后的圖像；H(S)代表目標圖像的熵，H(T*A)代表變換后圖譜的熵，H(S，T*A)代表目標圖像和變換后圖譜間的聯(lián)合熵；nMI是英文的歸一化互信息的縮寫(normalized mutual information)，nMI(S，T*A)表示S與T*A之間的歸一化互信息；空間變換T為包含9個參數(shù)的仿射變換T＝(p，q，r，u，v，w，φ，ω，θ)；其中，p，q，r分別為空間變換后x，y，z方向的位移；r，u，v分別為空間變換后x，y，z方向的比例縮放；φ，ω，θ分別為空間變換后繞x，y，z軸的轉角；若變換前點坐標為X＝(x，y，z)，則變換后的坐標為X′＝X′，y′，z′)，x′y′z′1=100p010q001r0001×cosφsinφ00-sinφcosφ0000100001×00-sinω00100sinω0cosω00001]]>×10000cosθsinθ00-sinθcosθ00001×u0000v0000w00001×xyz1]]>至此，定義了目標圖像S和經過空間變換T后的圖譜A之間的歸一化互信息；所述基于歸一化互信息的配準方法依次按以下步驟實現(xiàn)步驟(1.1).隨機確定初始的空間變換T0＝(T0，q0，r0，u0，v0，w0，φ0，ω0，θ0，)中九個參數(shù)值，并用由此得到的空間變換作用于圖譜A，得到變換后的圖譜T0*A；步驟(1.2).計算目標圖像S和經過空間變換T后的圖譜T*A間的聯(lián)合歸一化聯(lián)合直方圖h(l，k)h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}其中，S(x)＝l表示在目標圖像中像素x對應的灰度值為l，T*A(x)＝k表示在經過空間變換的圖譜T*A像素中x對應的灰度值為k；Num為像素x的數(shù)目；
步驟(1.3).計算目標圖像S和經過空間變換T后的圖譜T*A間的聯(lián)合概率分布pS，T*A(l，K)pS,T*A(l,k)=h(l,k)Σl,kh(l,k).]]>步驟(1.4).分別計算圖像S的邊緣概率分布pS(l)和T*A的邊緣概率分布pT*A(k)pS(l)=ΣkpS,T*A(l,k),]]>pT*A(k)=ΣlpS,T*A(l,k).]]>邊緣概率分布是概率統(tǒng)計理論中的概念，對于兩個隨機變量x，y組成的隨機向量(x，y)，其分量x的概率分布稱為隨機向量(x，y)的關于x的邊緣分布；這里是把圖像S和變換后的圖譜T*A分別作為兩個隨機變量l和k，然后分別計算l和k的邊緣概率分布；步驟(1.5).按以下公式分別計算目標圖像的熵H(S)，變換后圖譜的熵H(T*A)以及目標圖像和變換后圖譜間的聯(lián)合熵H(S，T*A)；H(S)=-ΣlpS(l)logpS(l)]]>H(T*A)=-ΣkpT*A(k)logpT*A(k)]]>H(S,T*A)=-ΣkpS,T*A(l,k)logS,T*A(l,k)]]>步驟(1.6).按下式計算目標圖像S和經過T變換后的圖譜T*A間的歸一化互信息nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>步驟(1.7).將歸一化互信息nMI(S，T*A)帶入Powell優(yōu)化算法對空間變換T進行反復優(yōu)化迭代，一直到收斂為止；Powell優(yōu)化算法是一種傳統(tǒng)的確定性優(yōu)化方法，并且被認為是目前最有效的無需求導數(shù)的優(yōu)方法；該算法的基本思想是，對于n維極值問題，首先沿著n個坐標方向求極小，經過n次之后得到n個共軛方向，然后沿n個共軛方向求極小，經過多次迭代后便可求得極小值；在這里，利用Powell優(yōu)化算法對空間變換T的值進行優(yōu)化，以求得到互信息nMI(S，T*A)的極小值；空間變換T含有九個分量，即(p，q，r，u，v，w，φ，ω，θ)，這九個分量分別對應Powell優(yōu)化算法中的九個維度，即利用Powell優(yōu)化算法求九維極值問題；當Powell優(yōu)化算法搜索結束，即算法收斂后，得到的nMI(S，T*A)的極小值所對應的空間變換值即為最優(yōu)空間變換T*，T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；則T**A為整體配準所要得到的配準圖譜；步驟(2).依次按以下步驟對感興趣的每一個器官進行單獨配準，確定每個感興趣的器官在目標圖像S中的大致對應區(qū)域；具體分為以下幾步步驟(2.1).按下式為每個感興趣器官確定灰度范圍區(qū)間RkRk＝[μk-λσk，μk+λσk]，其中，μk，σk為變換后圖譜T*A的第k個感興趣器官在目標圖像S中所占據(jù)的區(qū)域內的灰度均值和標準方差；λ為一個常數(shù)，用來調整灰度均勻范圍；試驗發(fā)現(xiàn)，λ的取值范圍為1.2～1.3；步驟(2.2).將步驟(1)中得到的最終的空間變換T*作對為每個器官單獨進行配準的初始空間變換T′，即令T′＝T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；步驟(2.3).按下式計算變換后的器官區(qū)域的相似性測度Mk(T′)Mk(T′)=Nkin(T′)-Nkout(T′).]]>其中，下標k代表對第k個器官進行配準，Nkin(T′)代表采用空間變換T′時器官k在配準圖譜中對應區(qū)域內包含的器官像素數(shù)目，Nkout(T′)代表該對應區(qū)域內包含的非器官像素數(shù)目；Nkin(T′)=N{x|x∈Ck(T′),S(x)∈Rk}Nkout(T′)=N{x|x∈Ck(T′),S(x)∉Rk}]]>其中，N{}表示計算集合中包含的像素數(shù)目，Ck(T′)表示器官k在空間變換T下對應的區(qū)域，S(x)表示像素x在CT或MRI圖像中的灰度值；步驟(2.4).將變換后的器官區(qū)域的相似性測度Mk(T′)帶入Powell優(yōu)化算法對T′進行反復迭代，直到收斂為止；具體過程即利用Powell優(yōu)化算法對空間變換T′的值進行優(yōu)化，以求得到相似性測度Mk(T′)的極小值；空間變換T′含有九個分量，即(p′，q′，r′，u′，v′，w′，φ′，ω′，θ′)，這九個分量分別對應Powell優(yōu)化算法中的九個維度，即利用Powell優(yōu)化算法求九維極值問題；當Powell優(yōu)化算法搜索結束，即算法收斂后，得到的Mk(T′)的極小值所對應的空間變換值即為最優(yōu)空間變換T′*，T′*＝(p′*，q′*，r′*，u′*，v′*，w′*，ω′*，θ′*)；則Ck(T′*)為單獨配準所要得到的器官區(qū)域；步驟(3).模糊連接分割；具體分為以下幾步步驟(3.1).設Ck(T′*)為模糊連接對第k個器官進行分割的初始區(qū)域，按下式計算該區(qū)域中的灰度直方圖Hk(l)Hk(l)＝N{x|x∈Ck，S(x)＝l}，l＝0，1……LHk(l)可理解為Ck(T′*)區(qū)域中灰度為l的像素x的個數(shù)；步驟(3.2).采用曲線擬合的方法對直方圖進行平滑和去噪，得到擬合后的直方圖 H^k(l)=Σj=1najexp[-(l-bjcj)2],]]>其中，n代表高斯分量的數(shù)目，即所要分割的器官的數(shù)目；例如我們對肝臟和脾臟兩種器官進行分割，則n＝2；aj，bj，cj依次代表第j個高斯分布的峰值、中心和寬度；步驟(3.3).擬合后的直方圖曲線可以看作器官在目標圖像中的灰度概率密度函數(shù)；選擇對應最大灰度概率的至少一個像素作為器官的初始種子點；步驟(3.4).考察每一個種子點相鄰區(qū)域的灰度概率之和；相鄰區(qū)域定義為與作為中心像素的種子點距離小于一個像素單位的距離的像素集合；所述灰度概率為所研究器官區(qū)域內具有該灰度的像素數(shù)目與區(qū)域內像素總數(shù)的比值；具有最大區(qū)域灰度概率總和的點將被最終確定為種子點；步驟(3.5).以由步驟(3.2)得到的直方圖 所對應的器官區(qū)域灰度分布的概率密度函數(shù)和最終種子點的位置為初始參數(shù)，用下式所述的用模糊連接算法計算整幅圖像中每一像素點與器官最終種子點的模糊連接強度cμ(p,q)=maxρ(p,q)∈P(p,q)[minz∈ρ(p,q)Hk^(l(z))]]]>其中，p表示最終的種子點，q表示整幅圖像中任意一點，cμ(p，q)表示p與q之間的模糊連接強度；ρ(p，q)表示p與q之間的任意一條連通路徑，P(p，q)表示p與q之間所有連通路徑的集合，z∈ρ(p，q)表示Z是路徑ρ(p，q)上的任意一像素點，l(z)為像素Z在目標圖像中的灰度值， 表示灰度值l(z)在直方圖 中所對應的概率密度值；模糊連接強度表征的該像素屬于器官區(qū)域的隸屬度，將每個像素點的隸屬度作為灰度所得到的圖像稱為模糊連接圖像；步驟(3.6).為每個器官的模糊連接圖像指定一個最佳閾值；所謂閾值，就是用來將屬于目標器官的像素和背景像素區(qū)分開的一個數(shù)值，規(guī)定灰度大于或等于此閾值的像素為屬于目標器官的像素，灰度小于此閾值的像素為背景像素；所謂最佳閾值，就是由此閾值得到的器官像素所集合成的區(qū)域最符合真實的器官區(qū)域；最佳閾值是將模糊連接圖像轉化為最終的分割結果的重要參量；確定最佳閾值具體步驟如下首先將閾值設定為1，然后逐步降低閾值；用于每一次降低得到的閾值來分割模糊連接圖像，得到器官的分割區(qū)域；觀察分割區(qū)域的形狀隨閾值的改變，如果在某一次閾值降低時，形狀劇烈改變，則停止降低閾值，將上一閾值定為最佳閾值；關于分割區(qū)域的形狀改變的衡量辦法和最佳閾值的確定方法，進一步詳細解釋如下a).分割區(qū)域的形狀改變的衡量的辦法采用面積相對變化值和緊密度相對變化值對來衡量分割區(qū)域的形狀改變；緊密度(Compactness)的定義為 面積相對變化值和緊密度相對變化值定義如下對于每次閾值降低，分別計算降低前后的器官面積和緊密度；令降低前的器官面積和緊密度為原面積和原緊密度，降低后的器官面積和緊密度為新面積和新緊密度，則面積和緊密度的相對變化定義為 b).用面積相對變化值和緊密度相對變化值確定最優(yōu)閾值的方法在閾值降低過程中，如果在某一次降低時，發(fā)現(xiàn)面積和緊密度的相對變化值大于預先設定的參考值ρ，則認為形狀發(fā)生了劇烈改變，從而停止閾值搜索，并將降低前的閾值定為最佳閾值；試驗表明，參考值ρ取值0.15～0.2能夠取得較準確的最佳閾值；步驟(3.7)、按照下式用步驟(3.6)得到的最佳閾值對模糊連接圖像進行分割，得到器官的二值分割圖像設第k個器官的最佳閾值為Tk，則該器官的二值分割結果為 其中，x代表任一像素，F(xiàn)k(x)為第k個器官的模糊連接圖像在像素x處的灰度值，Bk(x)為分割后的二值圖像在像素x處的灰度值；步驟(4).器官形狀修正，以解決因為人體腹腔內相互緊貼的兩個器官因為灰度相近而造成模糊連接分割可能產生的錯誤分割結果；醫(yī)學圖像中，經常出現(xiàn)兩個在人體腹腔中相互緊貼的器官在圖像中的灰度也相似的現(xiàn)象；對于這樣的兩個器官，醫(yī)學圖像中在其表面相貼處可能沒有明顯的邊界，這會使模糊連接算法認為這兩個器官是長在一起的，從而錯誤地分割成同一個器官；本研究提出以器官形狀穩(wěn)定為依據(jù)來檢測錯誤分割，并采用距離變換結合分水嶺分割算法的對錯誤分割結果進行形狀修正；具體流程如下步驟(4.1).對于多層醫(yī)學圖像，按以下步驟從上到下，采用逐層分析的方法檢測存在錯誤分割的層面；步驟(4.1.1).
設定當前層為算法正在分析的層面，當前層的序號為i下方相鄰層為與當前層相鄰的下方層面，下方相鄰層的序號為i-1；Oki為當前層中的器官分割區(qū)域；
Oki-1為下方相鄰層中器官分割區(qū)域；在當前層中，令器官分割區(qū)域內的各點到區(qū)域邊界的距離為負值，對于整幅圖像中所有不屬于該區(qū)域的點，即所有該區(qū)域以外的點，令其到區(qū)域邊界的距離為正值；步驟(4.1.2).對當前層i中的器官分割區(qū)域Oki進行距離變換，即對整幅圖像中每點按下式計算它們到區(qū)域邊界的最近距離Dki 其中，y代表了區(qū)域邊界的像素點中與x距離最近一點；Bki(x)為步驟(3.7)中得到的最終分割結果Bk*(x)在當前層i中的像素值；將距離值dki(x)作為灰度值賦給每一像素點，就得到了距離圖像Dki；步驟(4.1.3).提取下方相鄰層中區(qū)域Oki-1的邊界上各點像素的位置，則在Dki中的相同位置處的各點的Dki(x)值可認為是Oki-1邊界上的像素的距離值；得到所有Oki-1的邊界像素的距離值，計算這些距離值的標準方差，以確定器官在第i層和第i-1層分割結果間的相似程度；如果標準方差大于預先設定的參考值，則認為檢測到了分割錯誤，下一層即為存在錯誤分割的層面；根據(jù)經驗，這個參考值一般設定在5個像素距離；按照上述方法，從上到下檢測所有層面，得到所有存在錯誤分割的層面；步驟(4.2).按下式對存在錯誤分割的各層面的二值圖像進行距離變換；設定對于器官區(qū)域內所有像素點，令其距離值為負值；對于整幅圖像中區(qū)域外所有不屬于該區(qū)域的區(qū)域以外的點，令其距離值為負的極小值，即負無窮-∞；則距離圖Dki′(x) 其中y代表了與x距離最近的邊界點；Bki即存在錯誤分割的層面的二值圖像；步驟(4.3).按以下方法對步驟(4.2)得到的距離圖像Dki′應用分水嶺算法；步驟(4.3.1).對整幅圖像設定一個灰度閾值T；令T的初值為圖像Dki′中除了負無窮以外的最小值，即
T0=min({Dki(x)|Dki(x)≠-∞}),]]>其中，x代表圖像中任一個像素點的坐標；步驟(4.3.2).以d為步長逐漸增大閾值T，即Tn＝T0+n·d，其中，Tn表示第n次增加步長時閾值T的取值；n取值范圍1≤n≤K，其中K=(max(Dki(x))-T0)/d,]]>max(Dki(x))表示圖像Dki′的最大值；d取值固定，根據(jù)經驗取d＝3；步驟(4.3.3).利用第n次增加步長得到的閾值Tn按下式將圖像Dki′分割成二值圖像Wn 其中，x代表圖像中任一個像素點的坐標；Wn(x)＝1表示淹沒在水面以下的部分，Wn(x)＝0表示露在水面以上的陸地部分；對于前一次(即第n-1次)增加步長得到閾值Tn-1，按下式將圖像Dki′分割成二值圖像Wn-1 其中，x代表圖像中任一個像素點的坐標；步驟(4.3.4).在Tn的增長過程中， 每當Tn大于一個局部極大值，就說明圖像Wn-1內有兩個連通區(qū)域在Wn中合并成了一個連通域，則將Wn-1和Wn兩幅圖像進行異或操作，得到分水嶺圖像Sn 其中，x代表圖像中任一個像素點的坐標，表示異或操作，Sn(x)＝1表示像素x為分水嶺所在位置；步驟(4.3.5).對n取值范圍1≤n≤K中的每一個n值，均會得到一個相應的分水嶺圖像Sn；設所有的n值所對應的分水嶺圖像為S1，S2，…， 則分水嶺變換的結果圖像S為所有分水嶺圖像的集合
S=S1·S2·...·Smax(Dki(x));]]>圖像S中所有像素值為1的像素就表示圖像Dki′中的分水嶺所在的位置；由這些分水嶺即把圖像Dki′分割成若干個子區(qū)域；步驟(4.4).對分水嶺算法得到的圖像Dki′中的每個子區(qū)域，分別計算其與下方相鄰層面器官分割區(qū)域Oki-1的面積重合比率 把重合比率低于50％的區(qū)域視為錯誤分割區(qū)域，重新標記為背景區(qū)域，從而將其在距離圖像Dki′中除去。
實驗效果分析為了驗證本分割算法的可靠性，實驗選擇了5套來自于國內醫(yī)院的臨床CT腹部掃描圖像。這5套數(shù)據(jù)包括男性和女性，以及不同年齡的對象。從分割結果看出，試驗數(shù)據(jù)中的感興趣器官均得到了有效分割。配準后圖譜給出了大致分布區(qū)域，經過進一步的分割和修正，最終分割結果和手工分割結果相比，二者基本一致。


圖1為發(fā)明采用的VIP-Man圖譜中的一個層面。
圖2為CT圖像的一個層面。
圖3為計算機程序的主體流程圖。
圖4為將圖譜向CT圖像整體配準后，提取圖譜中感興趣器官的輪廓疊加顯示在CT圖像上。感興趣器官包括肝臟、腎臟和脾臟。
圖5為器官配準后，圖譜中感興趣器官的輪廓和CT圖像中器官位置對比。
圖6為各個器官的灰度直方圖和曲線擬合示意圖。
圖7為肝臟的模糊連接分割過程。從左至右依次為，CT圖像中肝臟局部，肝臟的模糊連接圖像，對模糊連接圖像進行閾值分割后的肝臟分割結果。
圖8為分水嶺算法示意圖。
圖9為對腎臟進行形狀修正過程。從左至右依次為，CT圖像中腎臟局部，腎臟的模糊連接分割結果，腎臟經過分水嶺算法被分割為多個區(qū)域和白色輪廓表示的相鄰層面腎臟區(qū)域對比，消除重合比率過低的區(qū)域后得到的腎臟修正后分割結果。
圖10對CT圖像的整體分割結果。
具體實施例方式
本發(fā)明提出了一種基于圖譜匹配的腹部器官分割方法，實現(xiàn)醫(yī)學圖像中腹部多器官的自動化分割。發(fā)明采用由美國虛擬人CT數(shù)據(jù)經過手工分割建立的圖譜，對腹部器官進行分割。關于該圖譜建立方法的詳細說明可以參見Xu等人2000年在《Health Physics》雜志發(fā)表的論文《VIP-manAn image-based whole-body adult male model constructed from color photographs ofthe visible human project for multi-particle Monte Carlo calculations》。該圖譜的示例可見圖1。通過配準工作，得到的是一個和目標圖像大致對齊的圖譜，目標圖像即需要做器官分割的臨床圖像，其示例見圖2。配準的目的就是將圖譜通過形變和目標圖像大致對齊，這種對齊集中于感興趣器官的對齊，而不是圖像所有區(qū)域的對齊，這樣就大大減輕了配準工作的壓力和困難程度。從而允許在配準中選擇較簡單的空間變換，只包含較少的自由度。優(yōu)化過程也相應的提高了速度，增強了穩(wěn)定性。在分割算法上，選擇模糊連接分割算法，它具有較強的穩(wěn)定性，能夠適應復雜圖像的分割。模糊連接分割面臨的一個主要問題是各項必要參數(shù)的指定，包括初始的種子點，感興趣器官的灰度分布特征等。利用配準后圖譜，本發(fā)明提出了這些參數(shù)的自動化指定方法，避免了常用的人工或經驗指定。已發(fā)現(xiàn)，由于相鄰器官間的灰度重疊或重合現(xiàn)象，模糊連接分割的結果可能包含著錯誤的區(qū)域。為了消除這些錯誤區(qū)域，提高分割結果的準確率，本研究基于器官在上下兩張相鄰層面的CT圖像中具有形狀一致性的原理，提出了器官形狀修正方法，對分割結果做進一步的完善。
本方法的計算機程序的主體流程圖可參加圖3。
本發(fā)明方法通過下述步驟實現(xiàn)1.整體配準整體配準的目標是消除圖譜和目標圖像間的整體差異，包括成像位置以及身體形態(tài)(高矮胖瘦等)的差異。整體配準，尋找的是一種空間變換，使圖譜經過這一空間變換后，能夠在整體上與目標圖像大致對齊。整體配準中選擇相似變換(包含相對于x、y、z軸3個坐標軸的平移、旋轉和比例變換，共9個參數(shù))，以歸一化的互信息作為目標函數(shù)，采用Powell算法進行優(yōu)化。
圖譜和目標圖像間歸一化互信息(Normalized Mutual Information，以下縮寫為nMI)定義為nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>其中A代表圖譜，T*A代表對圖譜做空間變換T后的圖像，H(S)代表目標圖像的熵，H(T*A)代表變換后圖譜的熵，H(S，T*A)代表目標圖像和變換后圖譜間的聯(lián)合熵。
在空間變換上，選擇包含9個參數(shù)的仿射變換T＝(p，q，r，u，v，w，φ，ω，θ)。其中，p，q，r分別為空間變換后x，y，z方向的位移；u，v，w分別為空間變換后x，y，z方向的比例縮放；φ，ω，θ分別為空間變換后饒x，y，z軸的轉角；令變換前點坐標為X＝(x，y，z)，變換后的坐標為X′＝(x′，y′，z′)，則x′y′z′1=100p010q001r0001×cosφsinφ00-sinφcosφ0000100001×00-sinω00100sinω0cosω00001]]>×10000cosθsinθ00-sinθcosθ00001×u0000v0000w00001×xyz1]]>H(S)=-ΣlpS(l)logpS(l)]]>H(T*A)=-ΣkpT*A(k)logpT*A(k)]]>H(S,T*A)=-ΣkpS,T*A(l,k)logS,T*A(l,k)]]>這里，pS，T*A(l，k)是圖像S和T*A間的聯(lián)合概率分布，可以用歸一化的聯(lián)合直方圖表示pS,T*A(l,k)=h(l,k)Σl,kh(l,k).]]>ps(l)，pT*A(k)是圖像S和T*A間的邊緣概率分布，由以下公式計算pS(l)=ΣkpS,T*A(l,k),]]>pT*A(k)=ΣlpS,T*A(l,k).]]>邊緣概率分布是概率統(tǒng)計理論中的概念，對于兩個隨機變量x，y組成的隨機向量(x，y)，其分量x的概率分布稱為隨機向量(x，y)的關于x的邊緣分布；這里是把圖像S和變換后的圖譜T*A分別作為兩個隨機變量l和k，然后分別計算l和k的邊緣概率分布；歸一化聯(lián)合直方圖h(l，k)定義為
h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}，即，像素在圖像S中灰度為l，并且對應像素在圖像T*A中為k的像素數(shù)目統(tǒng)計。
至此，定義了目標圖像S和經過空間變換T后的圖譜A之間的互信息。
Powell優(yōu)化算法是一種傳統(tǒng)的確定性優(yōu)化方法，并且被認為是目前最有效的無需求導數(shù)的優(yōu)方法。該算法的基本思想是，對于n維極值問題，首先沿著n個坐標方向求極小，經過n次之后得到n個共軛方向，然后沿n個共軛方向求極小，經過多次迭代后便可求得極小值。
在這里，利用Powell優(yōu)化算法對空間變換T的值進行優(yōu)化，以求得到互信息nMI(S，T*A)的極小值。空間變換T含有九個分量，即(p，q，r，u，v，w，φ，ω，θ)，這九個分量分別對應Powell優(yōu)化算法中的九個維度，即利用Powell優(yōu)化算法求九維極值問題。當Powell優(yōu)化算法搜索結束，即算法收斂后，得到的nMI(S，T*A)的極小值所對應的空間變換值即為最優(yōu)空間變換T*，T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；則T**A為整體配準所要得到的配準圖譜；圖4展示了將圖譜向CT圖像整體配準后的效果。我們將經過最優(yōu)空間T*變換后的圖譜中的感興趣器官的輪廓提取出來并疊加顯示在CT圖像上，以示效果。
2.器官配準器官配準的基本思想是進一步消除整體配準中未能解決的對應器官間的位置差異?；痉椒ㄊ菍D譜中每個感興趣器官分別地向個體圖像上配準，同時還要保證各個器官間不發(fā)生位置沖突。
對各個器官的單獨配準仍然選擇相似變換，采用Powell優(yōu)化算法。器官配準定義一個全新的目標函數(shù)，圖譜中的每一個感興趣器官在目標圖像中都存在一個對應區(qū)域，器官配準將努力使得該對應區(qū)域內包含盡可能多的器官像素，和盡可能少的非器官像素。由于這里還無法確定每個像素是否屬于某個器官，因此，采用一種推測標準。設共有n個感興趣器官，則對其中第k個感興趣器官，其灰度區(qū)間推測為，Rk＝[μk-λσk，μk+λσk]，其中，μk，σk為圖譜的第k個感興趣器官經配準后所占據(jù)的區(qū)域內的灰度均值和標準方差。
其中λ為一個常數(shù)，用來調整灰度均勻范圍。試驗發(fā)現(xiàn)，λ的取值范圍為1.2～1.3。如像素灰度在區(qū)間Rk內，則認為其屬于器官k像素；否則，認為其屬于非器官像素。配準的目標是通過最優(yōu)化算法，對每一個器官分別確定一個新的相似變換；這里采用的最優(yōu)化算法為Powell算法。新的變換使得器官對應區(qū)域內器官像素盡可能多，非器官像素盡可能少。定義器官配準的相似性測度Mk(T)如下Mk(T)=Nkin(T)-Nkout(T).]]>其中，下標k代表對第k個器官進行配準，Nkin(T)代表采用空間變換T時器官k的對應區(qū)域內包含的器官像素數(shù)目，Nkout(T)代表對應區(qū)域內包含的非器官像素數(shù)目。
Nkin(T)=N{x|x∈Ck(T),S(x)∈Rk}Nkout(T)=N{x|x∈Ck(T),S(x)∉Rk}]]>其中，N{}表示計算集合中包含的像素數(shù)目，Ck(T)表示器官k在空間變換T下對應的初始區(qū)域，S(x)表示像素在CT或MRI圖像中的灰度值。
Powell優(yōu)化算法的具體優(yōu)化過程為將變換后的器官區(qū)域的相似性測度Mk(T′)帶入Powell優(yōu)化算法，利用Powell優(yōu)化算法對空間變換T′的值進行優(yōu)化，以求得到相似性測度Mk(T′)的極小值；空間變換T′含有九個分量，即(p′，q′，r′，u′，v′，w′，φ′，ω′，θ′)，這九個分量分別對應Powell優(yōu)化算法中的九個維度，即利用Powell優(yōu)化算法求九維極值問題；當Powell優(yōu)化算法搜索結束，即算法收斂后，得到的Mk(T′)的極小值所對應的空間變換值即為最優(yōu)空間變換T′*，T′*＝(p′*，q′*，r′*，u′*，v′*，w′*，φ′*，ω′*，θ′*)；則Ck(T′*)為單獨配準所要得到的器官區(qū)域；圖5展示了經過器官配準后，圖譜中的感興趣器官的輪廓和CT圖像中器官位置對比。
3.模糊連接分割模糊連接算法由J.Udupa提出，它通過定義圖像中每一像素點與器官種子點的模糊連接強度來得到每個像素屬于該器官的隸屬度。將每個像素點的隸屬度作為灰度的圖像即模糊連接圖像。為了將模糊連接方法應用于器官分割，首先需要指定每個感興趣器官的灰度分布特征和種子點位置(器官的種子點是指模糊連接算法分割開始分割時的初始像素點)。這些參數(shù)確定之后，模糊連接算法會分別計算每個感興趣器官的模糊連接圖像；模糊連接圖像是灰度圖像，其中每個像素點的灰度值代表了該點與種子點的模糊連接程度。模糊連接圖像生成后，還需要進一步指定一個閾值，將其轉化為對器官的二值分割圖像，作為最終的分割結果。
配準后的圖譜給出了器官在目標圖像中的大致位置和范圍，稱其為器官的初始區(qū)域。統(tǒng)計初始區(qū)域中像素的灰度直方圖。以第k個器官為例，灰度直方圖的計算公式為，Hk(l)＝N{x | x∈Ck，S(x)＝l}， l＝0，1……L，即在器官k的初始區(qū)域中灰度為l的像素的數(shù)目。
器官灰度直方圖反映了像素灰度值與像素屬于器官的可能性大小之間的關系。由于初始區(qū)域中包含的不全部是器官像素，還有非器官像素。因此采用曲線擬合的方法對直方圖進行平滑和去噪。曲線擬合將直方圖視為若干個高斯函數(shù)的疊加，采用最小方差進行擬合。
器官擬合公式為H^k(l)=Σj=1najexp[-(l-bjcj)2],]]>其中，n代表高斯分量的數(shù)目，即所要分割的器官的數(shù)目；例如我們對肝臟和脾臟兩種器官進行分割，則n＝2；aj，bj，cj依次代表第j個高斯分布的峰值、中心和寬度。在高斯分量數(shù)目和擬合精度之間，需要進行平衡，用盡可能少的高斯分量取得盡可能高的擬合精度。這里，采用調整的殘值平方(Adjusted R-square，簡稱AR)來尋找最優(yōu)的高斯分量數(shù)目。
AR=1-SSE*LSST*(L-n*3),]]>其中SSE=Σl=0L(Hk(l)-H^k(l))2]]>SST=Σl=0L(Hk(l)-H‾k)2,]]>H‾ki=Σl=0LHk(l)/(L+1).]]>具體的高斯擬合過程為(1)令n＝0；(2)令n＝n+1；(3)采用最小SSE方法，擬合式 (4)計算AR；(5)當AR＜AR并且n＜5返回執(zhí)行步驟(2)；
在上述步驟中，TAR是依據(jù)經驗設定的閾值，用于控制平衡擬合精度和高斯分量數(shù)目。試驗表明，TAR＝0.95～0.97是一個比較合理的閾值設定。
圖6展示了各個器官的灰度直方圖和曲線擬合示意圖。擬合后的直方圖曲線可以看作器官在目標圖像中的灰度概率密度函數(shù)。選擇對應最大灰度概率的至少一個像素作為器官的初始種子點。對應最大灰度概率的像素數(shù)目可能大于1，再進一步考察像素相鄰區(qū)域的灰度概率之和。相鄰區(qū)域一般定義為與中心像素距離小于一個像素單位的像素集合。所述區(qū)域灰度概率為所研究器官區(qū)域內具有該灰度的像素數(shù)目與區(qū)域內像素總數(shù)的比值；具有最大區(qū)域灰度概率總和的點將被最終確定為種子點。同時規(guī)定，種子點必須位于器官的初始區(qū)域內部。
確定了器官的灰度分布特征和初始種子點，可以計算出每一個器官各自的模糊連接圖像。以上面得到的直方圖 所對應的器官區(qū)域灰度分布的概率密度函數(shù)和最終種子點的位置為初始參數(shù)，用下式所述的用模糊連接算法計算整幅圖像中每一像素點與器官最終種子點的模糊連接強度cμ(p,q)=maxρ(p,q)∈P(p,q)[minz∈ρ(p,q)H^k(l(z))]]]>其中p表示最終的種子點，q表示整幅圖像中任意一點，cμ(p，q)表示p與q之間的模糊連接強度；ρ(p，q)表示p與q之間的任意一條連通路徑，P(p，q)表示p與q之間所有連通路徑的集合，z∈ρ(p，q)表示Z是路徑ρ(p，q)上的任意一像素點，l(z)為像素Z在目標圖像中的灰度值， 表示灰度值l(z)在直方圖 中所對應的概率密度值；模糊連接強度表征的該像素屬于器官區(qū)域的隸屬度，將每個像素點的隸屬度作為灰度所得到的圖像稱為模糊連接圖像。而要將器官真正分割出來，還需要指定一個模糊連接強度閾值，對模糊連接圖像進行二值分割。所謂閾值，就是用來將目標像素和背景像素區(qū)分開的一個數(shù)值；在模糊連接圖像中，我們認為灰度值大于或等于閾值的像素為器官像素，灰度值小于閾值像素為背景像素。所謂最佳閾值，就是由此閾值得到的器官像素所集合成的區(qū)域最符合真實的器官區(qū)域；最佳閾值是將模糊連接圖像轉化為最終的分割結果的重要參量。
為了找到最佳的模糊連接分割閾值，提出如下的搜索方法。首先將閾值設定為1，然后逐步降低閾值；用于每一次降低得到的閾值來分割模糊連接圖像，得到器官的分割區(qū)域；分割方法如下設第k個器官的閾值為Tk，則該器官的閾值分割結果為
其中x代表任一像素，F(xiàn)k(x)為第k個器官的模糊連接圖像在像素x處的灰度值。
對于每一次閾值降低，觀察分割結果隨閾值的形狀改變。如果在某一次閾值降低時，形狀劇烈改變，則停止閾值搜索，將上一閾值定為最佳閾值。這里對形狀改變的衡量使用了兩個參量，即面積相對變化值和緊密度相對變化值。緊密度(Compactness)的定義為 對于每次閾值降低，分別計算降低前后的器官面積和緊密度；令降低前的器官面積和緊密度稱為原面積和原緊密度，降低后的器官面積和緊密度稱為新面積和新緊密度，則面積和緊密度的相對變化定義為 如果在某一次閾值降低時，發(fā)現(xiàn)面積和緊密度的相對變化值大于預先設定的參考值ρ，則認為形狀發(fā)生了劇烈改變，從而停止閾值搜索，并將降低前的閾值定為最佳閾值。試驗表明，參考值ρ取值0.15～0.2能夠取得較準確的最佳閾值。
接下來，用上面確定的最佳閾值將模糊連接圖像分割為二值圖像，作為最終的輸出。具體方法如下設第k個器官的最佳閾值為Tk*，則該器官的二值分割結果為 其中x代表任一像素，F(xiàn)k(x)為第k個器官的模糊連接圖像在像素x處的灰度值，Bk*(x)為分割后的二值圖像在像素x處的灰度值。
圖7以肝臟為例展示了模糊連接分割過程。從左至右依次為，CT圖像中肝臟局部，肝臟的模糊連接圖像，對模糊連接圖像進行閾值分割后的肝臟分割結果。
4.器官形狀修正醫(yī)學圖像中，經常出現(xiàn)兩個在人體腹腔中相互緊貼的器官在圖像中的灰度也相似的現(xiàn)象；對于這樣的兩個器官，醫(yī)學圖像中在其表面相貼處可能沒有明顯的邊界，這會使模糊連接算法認為這兩個器官是長在一起的，從而錯誤地分割成同一個器官；本研究提出以器官形狀穩(wěn)定為依據(jù)來檢測錯誤分割，并采用距離變換結合分水嶺分割算法的對錯誤分割結果進行形狀修正；具體流程如下
1).對于多層醫(yī)學圖像，采用逐層分析的方法檢測存在錯誤分割的層面；檢測順序從最上方層面到最下方層面；對于每次層面移動，稱算法正在分析的層面為當前層，與當前層相鄰的下方層面稱為下方相鄰層；設當前層的序號為i，則下方相鄰層的序號為i-1；在當前層中，設Oki代表器官分割區(qū)域；對Oki進行距離變換，即計算區(qū)域內每個離散點到區(qū)域邊界的距離；對于區(qū)域內所有點，令其距離值為負值；對于區(qū)域外所有點(范圍包括整幅圖像中所有不屬于該區(qū)域的點)，令其距離值為正值，這樣得到距離圖像為Dki， 其中x代表圖像中任意一處的離散像素點，y代表了區(qū)域邊界上與x距離最近的點；Bki(x)為模糊連接最終分割結果Bk*(x)在第當前層的像素值。
在下方相鄰層中，假定Oki-1代表器官分割區(qū)域；提取Oki-1的邊界像素的位置，則在Dki中的相同位置處的Dki(x)值可認為是Oki-1邊界上的像素與Oki邊界上的對應像素的距離值；當?shù)玫剿蠴ki-1的邊界像素的這種距離值后，計算這些距離值的標準方差，該標準方差反映了器官在第i層和第i-1層分割結果間的相似程度。如果器官從一層到另一層只是均勻的膨脹或收縮，則標準方差取值應小于預定義參考值。因此，可以采用器官在相鄰層面間邊界距離標準方差的大小來衡量是否出現(xiàn)了錯誤的異常突出區(qū)域。如果標準方差大于參考值，則認為檢測到了分割錯誤。根據(jù)試驗經驗，這個參考值一般設定在5個像素距離。
按照上述方法，從上到下檢測所有層面，得到所有存在錯誤分割的層面。
2).對存在錯誤分割的層面的二值圖像進行距離變換；這里的距離變換與上面有所不同；對于區(qū)域內所有點，令其距離值為負值；對于區(qū)域外所有點(范圍包括整幅圖像中所有不屬于該區(qū)域的點)，令其距離值負的極小值；這樣得到距離圖像Dki′為 其中x代表圖像中任意一處的離散像素點，y代表了區(qū)域邊界上與x距離最近的點；Bki即存在錯誤分割的層面的二值圖像；3).對圖像Dki′應用分水嶺算法。
借助漲水法來介紹分水嶺技術。分水嶺(Watershed)算法把圖像中像素的灰度值看作地形圖上的海拔高度，則圖像中灰度高的部分對應山峰，灰度低的部分對應山谷。利用圖8(a)來進行討論，為了簡便，僅畫出圖像的一維剖面圖。橫坐標代表圖像中的像素點，縱坐標代表像素的灰度值。假設有水從谷底涌出并水位逐漸升高。如果從兩個谷底涌出的水水位高過了其間的山峰，這些水就會匯合。如果要阻止水的匯合，則需要在山峰上修建水壩，且水壩的高度隨水位的上升而增高。這個過程隨著全部山峰都被水淹沒而結束。修建的所有水壩把整幅地形圖分割成很多區(qū)域，這些水壩就構成了這片土地的分水嶺。被分水嶺分割出來的每一個區(qū)域對應于圖像中的一個子區(qū)域，因此分水嶺算法將距離圖像Dki′分成多個子區(qū)域。用分水嶺算法將距離圖像Dki′分成多個子區(qū)域的具體過程如下仍以圖8(a)而言，圖中出現(xiàn)的三個山峰表示從圖像Dki′中任意取出的一部分，山峰表示灰度高的像素值，山谷表示灰度低的像素值。設灰度閾值T代表水面的高度；令T的初值為圖像Dki′中除負無窮以外的最小值，即T0=min({Dki(x)|Dki(x)≠-∞}),]]>其中，x代表圖像中任一個像素點的坐標。
以d為步長逐漸增大閾值T，即Tn＝T0+n·d，其中，Tn表示第n次增加步長時閾值T的取值；n取值范圍1≤n≤K；K=(max(Dki(x))-T0)/d,]]>max(Dki(x))表示圖像Dki′的最大值；d取值固定，根據(jù)經驗取d＝3。
利用第n次增加步長得到的閾值Tn(如圖8(c)所示)，可以按下式將圖像Dki′分割成二值圖像Wn 其中，x代表圖像中任一個像素點的坐標。Wn(x)＝1表示淹沒在水面以下的部分，Wn(x)＝0表示露在水面以上的陸地部分。
對于前一次(即第n-1次)增加步長得到閾值Tn-1(如圖8(b)所示)，可以按下式將圖像Dki′分割成二值圖像Wn-1
其中，x代表圖像中任一個像素點的坐標。
在Tn的增長過程中，每當Tn大于一個局部極大值，就說明圖像Wn-1內的兩個連通區(qū)域(如圖8(b)中的連通域Cj和Cj+1)在Wn中合并成了一個連通域(如圖8(c)中的Cj，j+1)，則將Wn-1和Wn兩幅圖像進行異或操作，得到分水嶺圖像Sn 其中，x代表圖像中任一個像素點的坐標，表示異或操作，Sn(x)＝1表示像素x為分水嶺所在位置，如圖8(d)所示。
因為n取值范圍1≤n≤ K，對于這個范圍的每一個n值，均會得到一個相應的分水嶺圖像Sn；則所有的n值所對應的分水嶺圖像為S1，S2，…， 分水嶺變換的結果圖像S為所有分水嶺圖像的集合S=S1·S2·...·Smax(Dki(x)).]]>則圖像S中所有像素值為1的像素就表示圖像Dki′中的分水嶺所在的位置；由這些分水嶺即把圖像Dki′分割成若干個子區(qū)域。
4).對分水嶺算法得到的每個子區(qū)域，分別計算其與下方相鄰層面器官分割區(qū)域Oki-1的面積重合比率；面積重合比率定義為 基于器官在相鄰層面中形狀一致的原則，認為重合比率低于50％的區(qū)域為應為錯誤分割區(qū)域，應重新標記為背景區(qū)域，，從而將其在距離圖像Dki′中除去。
圖9以腎臟為例說明了形狀修正過程。從左至右依次為，CT圖像中腎臟局部，腎臟的模糊連接分割結果，腎臟經過分水嶺算法被分割為多個區(qū)域和白色輪廓表示的下方相鄰層面腎臟區(qū)域對比，消除重合比率過低的區(qū)域后得到的腎臟修正后分割結果。
至此，就全部完成了對腹部醫(yī)學圖像的分割，最終的分割結果可見圖10。
權利要求
1.一種醫(yī)學圖像中腹部器官分割的方法，具體特征在于，該方法是在PC機中依次按以下步驟實現(xiàn)的算法的主要步驟包括步驟(1)，在PC機中，采用基于歸一化互信息的配準方法，將圖譜向通過CT或核磁掃描得到的個體圖像作整體配準，以消除圖譜和目標圖像間的整體差異；所述圖譜是指對已有的醫(yī)學圖象經解剖學專家手工分割后對每個器官區(qū)域做好標記的圖象，它可以為PC機提供人體解剖結構的參考信息；所述的歸一化互信息nMI(S，T*A)用下式表示nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>其中，S為目標圖像，T*A代表對圖譜A做空間變換T后的圖像；H(S)代表目標圖像的熵，H(T*A)代表變換后圖譜的熵，H(S，T*A)代表目標圖像和變換后圖譜間的聯(lián)合熵；nMI是英文的歸一化互信息的縮寫(normalized mutual information)，nMI(S，T*A)表示S與T*A之間的歸一化互信息；空間變換T為包含9個參數(shù)的仿射變換T＝(p，q，r，u，v，w，φ，ω，θ)；其中，p，q，r分別為空間變換后x，y，z方向的位移；r，u，v分別為空間變換后x，y，z方向的比例縮放；φ，ω，θ分別為空間變換后繞x，y，z軸的轉角；若變換前點坐標為X＝(x，y，z)，則變換后的坐標為X′＝(x′，y′，z′)，x′y′z′1=100p010q001r0001×cosφsinφ00-sinφcosφ0000100001×00-sinω00100sinω0cosω00001]]>×10000cosθsinθ00-sinθcosθ00001×u0000v0000w00001×xyz1]]>至此，定義了目標圖像S和經過空間變換T后的圖譜A之間的歸一化互信息；所述基于歸一化互信息的配準方法依次按以下步驟實現(xiàn)步驟(1.1).隨機確定初始的空間變換T0＝(p0，q0，r0，u0，v0，w0，φ0，ω0，θ0)中九個參數(shù)值，并用由此得到的空間變換作用于圖譜A，得到變換后的圖譜T0*A；步驟(1.2).計算目標圖像S和經過空間變換T后的圖譜T*A間的聯(lián)合歸一化聯(lián)合直方圖h(l，k)，h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}其中，S(x)＝l表示在目標圖像中像素x對應的灰度值為l，T*A(x)＝k表示在經過空間變換的圖譜T*A像素中x對應的灰度值為k；Num為像素x的數(shù)目；步驟(1.3).計算目標圖像S和經過空間變換T后的圖譜T*A間的聯(lián)合概率分布pS，T*A(l，k)pS,T*A(l,k)=h(l,k)Σl,kh(l,k).]]>步驟(1.4).分別計算圖像S的邊緣概率分布pS(l)和T*A的邊緣概率分布pT*A(k)PS(l)=ΣkPS,T*A(l,k),]]>pT*A(k)=ΣlpS,T*A(l,k).]]>邊緣概率分布是概率統(tǒng)計理論中的概念，對于兩個隨機變量x，y組成的隨機向量(x，y)，其分量x的概率分布稱為隨機向量(x，y)的關于x的邊緣分布；這里是把圖像S和變換后的圖譜T*A分別作為兩個隨機變量l和k，然后分別計算l和k的邊緣概率分布；步驟(1.5).按以下公式分別計算目標圖像的熵H(S)，變換后圖譜的熵H(T*A)以及目標圖像和變換后圖譜間的聯(lián)合熵H(S，T*A)；H(S)=-ΣlpS(l)logpS(l)]]>H(T*A)=-ΣkpT*A(k)logpT*A(k)]]>H(S,T*A)=-ΣkpS,T*A(l,k)logS,T*A(l,k)]]>步驟(1.6).按下式計算目標圖像S和經過T變換后的圖譜T*A間的歸一化互信息nMI(S,T*A)=H(S)+H(T*A)H(S,T*A)]]>步驟(1.7).將歸一化互信息nMI(S，T*A)帶入Powell優(yōu)化算法對空間變換T進行反復優(yōu)化迭代，一直到收斂為止；Powell優(yōu)化算法是一種傳統(tǒng)的確定性優(yōu)化方法，并且被認為是目前最有效的無需求導數(shù)的優(yōu)方法；該算法的基本思想是，對于n維極值問題，首先沿著n個坐標方向求極小，經過n次之后得到n個共軛方向，然后沿n個共軛方向求極小，經過多次迭代后便可求得極小值；在這里，利用Powell優(yōu)化算法對空間變換T的值進行優(yōu)化，以求得到互信息nMI(S，T*A)的極小值；空間變換T含有九個分量，即(p，q，r，u，v，w，φ，ω，θ)，這九個分量分別對應Powell優(yōu)化算法中的九個維度，即利用Powell優(yōu)化算法求九維極值問題；當Powell優(yōu)化算法搜索結束，即算法收斂后，得到的nMI(S，T*A)的極小值所對應的空間變換值即為最優(yōu)空間變換T*，T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；則T**A為整體配準所要得到的配準圖譜；步驟(2).依次按以下步驟對感興趣的每一個器官進行單獨配準，確定每個感興趣的器官在目標圖像S中的大致對應區(qū)域；具體分為以下幾步步驟(2.1).按下式為每個感興趣器官確定灰度范圍區(qū)間RkRk＝[μk-λσk，μk+λσk]，其中，μk，σk為變換后圖譜T*A的第k個感興趣器官在目標圖像S中所占據(jù)的區(qū)域內的灰度均值和標準方差；λ為一個常數(shù)，用來調整灰度均勻范圍；試驗發(fā)現(xiàn)，λ的取值范圍為1.2～1.3；步驟(2.2).將步驟(1)中得到的最終的空間變換T*作對為每個器官單獨進行配準的初始空間變換T′，即令T′＝T*＝(p*，q*，r*，u*，v*，w*，φ*，ω*，θ*)；步驟(2.3).按下式計算變換后的器官區(qū)域的相似性測度Mk(T′)Mk(T′)=Nkin(T′)-Nkout(T′).]]>其中，下標k代表對第k個器官進行配準，Nkin(T′)代表采用空間變換T′時器官k在配準圖譜中對應區(qū)域內包含的器官像素數(shù)目，Nkout(T′)代表該對應區(qū)域內包含的非器官像素數(shù)目；Nkin(T′)=N{x|x∈Ck(T′),S(x)∈Rk}Nkout(T′)=N{x|x∈Ck(T′),S(x)∉Rk}]]>其中，N{}表示計算集合中包含的像素數(shù)目，Ck(T′)表示器官k在空間變換T下對應的區(qū)域，S(x)表示像素x在CT或MRI圖像中的灰度值；步驟(2.4).將變換后的器官區(qū)域的相似性測度Mk(T′)帶入Powell優(yōu)化算法對T′進行反復迭代，直到收斂為止；具體過程即利用Powell優(yōu)化算法對空間變換T′的值進行優(yōu)化，以求得到相似性測度Mk(T′)的極小值；空間變換T′含有九個分量，即(p′，q′，r′，u′，v′，w′，φ′，ω′，θ′)，這九個分量分別對應Powell優(yōu)化算法中的九個維度，即利用Powell優(yōu)化算法求九維極值問題；當Powell優(yōu)化算法搜索結束，即算法收斂后，得到的Mk(T′)的極小值所對應的空間變換值即為最優(yōu)空間變換T*′，T′*＝(p′*，q′*，r′*，u′*，v′*，w′*，φ′*，ω′*，θ′*)；則Ck(T′*)為單獨配準所要得到的器官區(qū)域；步驟(3).模糊連接分割；具體分為以下幾步步驟(3.1).設Ck(T′*)為模糊連接對第k個器官進行分割的初始區(qū)域，按下式計算該區(qū)域中的灰度直方圖Hk(l)Hk(l)＝N{x|x∈Ck，S(x)＝l}，l＝0，1……LHk(l)可理解為Ck(T′*)區(qū)域中灰度為l的像素x的個數(shù)；步驟(3.2).采用曲線擬合的方法對直方圖進行平滑和去噪，得到擬合后的直方圖 H^k(l)=Σj=1najexp[-(l-bjcj)2],]]>其中，n代表高斯分量的數(shù)目，即所要分割的器官的數(shù)目；例如我們對肝臟和脾臟兩種器官進行分割，則n＝2；aj，bj，cj依次代表第j個高斯分布的峰值、中心和寬度；步驟(3.3).擬合后的直方圖曲線可以看作器官在目標圖像中的灰度概率密度函數(shù)；選擇對應最大灰度概率的至少一個像素作為器官的初始種子點；步驟(3.4).考察每一個種子點相鄰區(qū)域的灰度概率之和；相鄰區(qū)域定義為與作為中心像素的種子點距離小于一個像素單位的距離的像素集合；所述灰度概率為所研究器官區(qū)域內具有該灰度的像素數(shù)目與區(qū)域內像素總數(shù)的比值；具有最大區(qū)域灰度概率總和的點將被最終確定為種子點；步驟(3.5).以由步驟(3.2)得到的直方圖 所對應的器官區(qū)域灰度分布的概率密度函數(shù)和最終種子點的位置為初始參數(shù)，用下式所述的用模糊連接算法計算整幅圖像中每一像素點與器官最終種子點的模糊連接強度cμ(p,q)=maxρ(p,q)∈p(p,q)[minz∈ρ(p,q)H^k(l(z))]]]>其中，p表示最終的種子點，q表示整幅圖像中任意一點，cμ(p，q)表示p與q之間的模糊連接強度；ρ(p，q)表示p與q之間的任意一條連通路徑，P(p，q)表示p與q之間所有連通路徑的集合，z∈ρ(p，q)表示Z是路徑ρ(p，q)上的任意一像素點，l(z)為像素Z在目標圖像中的灰度值， 表示灰度值l(z)在直方圖 中所對應的概率密度值；模糊連接強度表征的該像素屬于器官區(qū)域的隸屬度，將每個像素點的隸屬度作為灰度所得到的圖像稱為模糊連接圖像；步驟(3.6).為每個器官的模糊連接圖像指定一個最佳閾值；所謂閾值，就是用來將屬于目標器官的像素和背景像素區(qū)分開的一個數(shù)值，規(guī)定灰度大于或等于此閾值的像素為屬于目標器官的像素，灰度小于此閾值的像素為背景像素；所謂最佳閾值，就是由此閾值得到的器官像素所集合成的區(qū)域最符合真實的器官區(qū)域；最佳閾值是將模糊連接圖像轉化為最終的分割結果的重要參量；確定最佳閾值具體步驟如下首先將閾值設定為1，然后逐步降低閾值；用于每一次降低得到的閾值來分割模糊連接圖像，得到器官的分割區(qū)域；觀察分割區(qū)域的形狀隨閾值的改變，如果在某一次閾值降低時，形狀劇烈改變，則停止降低閾值，將上一閾值定為最佳閾值；關于分割區(qū)域的形狀改變的衡量辦法和最佳閾值的確定方法，進一步詳細解釋如下a).分割區(qū)域的形狀改變的衡量的辦法采用面積相對變化值和緊密度相對變化值對來衡量分割區(qū)域的形狀改變；緊密度(Compactness)的定義為 面積相對變化值和緊密度相對變化值定義如下對于每次閾值降低，分別計算降低前后的器官面積和緊密度；令降低前的器官面積和緊密度為原面積和原緊密度，降低后的器官面積和緊密度為新面積和新緊密度，則面積和緊密度的相對變化定義為 b).用面積相對變化值和緊密度相對變化值確定最優(yōu)閾值的方法在閾值降低過程中，如果在某一次降低時，發(fā)現(xiàn)面積和緊密度的相對變化值大于預先設定的參考值ρ，則認為形狀發(fā)生了劇烈改變，從而停止閾值搜索，并將降低前的閾值定為最佳閾值；試驗表明，參考值ρ取值0.15～0.2能夠取得較準確的最佳閾值；步驟(3.7)、按照下式用步驟(3.6)得到的最佳閾值對模糊連接圖像進行分割，得到器官的二值分割圖像設第k個器官的最佳閾值為Tk，則該器官的二值分割結果為 其中，x代表任一像素，F(xiàn)k(x)為第k個器官的模糊連接圖像在像素x處的灰度值，Bk(x)為分割后的二值圖像在像素x處的灰度值；步驟(4).器官形狀修正，以解決因為人體腹腔內相互緊貼的兩個器官因為灰度相近而造成模糊連接分割可能產生的錯誤分割結果；醫(yī)學圖像中，經常出現(xiàn)兩個在人體腹腔中相互緊貼的器官在圖像中的灰度也相似的現(xiàn)象；對于這樣的兩個器官，醫(yī)學圖像中在其表面相貼處可能沒有明顯的邊界，這會使模糊連接算法認為這兩個器官是長在一起的，從而錯誤地分割成同一個器官；本研究提出以器官形狀穩(wěn)定為依據(jù)來檢測錯誤分割，并采用距離變換結合分水嶺分割算法的對錯誤分割結果進行形狀修正；具體流程如下步驟(4.1).對于多層醫(yī)學圖像，按以下步驟從上到下，采用逐層分析的方法檢測存在錯誤分割的層面；步驟(4.1.1).設定當前層為算法正在分析的層面，當前層的序號為i下方相鄰層為與當前層相鄰的下方層面，下方相鄰層的序號為i-1；Oki為當前層中的器官分割區(qū)域；Oki-1為下方相鄰層中器官分割區(qū)域；在當前層中，令器官分割區(qū)域內的各點到區(qū)域邊界的距離為負值，對于整幅圖像中所有不屬于該區(qū)域的點，即所有該區(qū)域以外的點，令其到區(qū)域邊界的距離為正值；步驟(4.1.2).對當前層i中的器官分割區(qū)域Oki進行距離變換，即對整幅圖像中每點按下式計算它們到區(qū)域邊界的最近距離Dki 其中，y代表了區(qū)域邊界的像素點中與x距離最近一點；Bik(x)為步驟(3.7)中得到的最終分割結果Bk*(x)在當前層i中的像素值；將距離值Dki(x)作為灰度值賦給每一像素點，就得到了距離圖像Dki；步驟(4.1.3).提取下方相鄰層中區(qū)域Oki-1的邊界上各點像素的位置，則在Dki中的相同位置處的各點的Dki(x)值可認為是Oki-1邊界上的像素的距離值；得到所有Oki-1的邊界像素的距離值，計算這些距離值的標準方差，以確定器官在第i層和第i-1層分割結果間的相似程度；如果標準方差大于預先設定的參考值，則認為檢測到了分割錯誤，下一層即為存在錯誤分割的層面；根據(jù)經驗，這個參考值一般設定在5個像素距離；按照上述方法，從上到下檢測所有層面，得到所有存在錯誤分割的層面；步驟(4.2).按下式對存在錯誤分割的各層面的二值圖像進行距離變換；設定對于器官區(qū)域內所有像素點，令其距離值為負值；對于整幅圖像中區(qū)域外所有不屬于該區(qū)域的區(qū)域以外的點，令其距離值為負的極小值，即負無窮-∞；則距離圖像Dki′(x) 其中y代表了與x距離最近的邊界點；Bki即存在錯誤分割的層面的二值圖像；步驟(4.3).按以下方法對步驟(4.2)得到的距離圖像Dki′應用分水嶺算法；步驟(4.3.1).對整幅圖像設定一個灰度閾值T；令T的初值為圖像Dki′中除了負無窮以外的最小值，即T0=min({Dki(x)|Dki(x)≠-∞}),]]>其中，x代表圖像中任一個像素點的坐標；步驟(4.3.2).以d為步長逐漸增大閾值T，即Tn＝T0+n·d，其中，Tn表示第n次增加步長時閾值T的取值； n取值范圍1≤n≤K，其中K=(max(Dki(x))-T0)/d]]>，max(Dki(x))表示圖像Dki′的最大值；d取值固定，根據(jù)經驗取d＝3；步驟(4.3.3).利用第n次增加步長得到的閾值Tn按下式將圖像Dki′分割成二值圖像Wn 其中，x代表圖像中任一個像素點的坐標；Wn(x)＝1表示淹沒在水面以下的部分，Wn(x)＝0表示露在水面以上的陸地部分；對于前一次(即第n-1次)增加步長得到閾值Tn-1，按下式將圖像Dki′分割成二值圖像Wn-1 其中，x代表圖像中任一個像素點的坐標；步驟(4.3.4).在Tn的增長過程中，每當Tn大于一個局部極大值，就說明圖像Wn-1內有兩個連通區(qū)域在Wn中合并成了一個連通域，則將Wn-1和Wn兩幅圖像進行異或操作，得到分水嶺圖像Sn 其中，x代表圖像中任一個像素點的坐標，表示異或操作，Sn(x)＝1表示像素x為分水嶺所在位置；步驟(4.3.5).對n取值范圍1≤n≤K中的每一個n值，均會得到一個相應的分水嶺圖像Sn；設所有的n值所對應的分水嶺圖像為S1，S2，…，Smax(Dki(x))，則分水嶺變換的結果圖像S為所有分水嶺圖像的集合S=S1·S2·...·Smax(Dki(x));]]>圖像S中所有像素值為1的像素就表示圖像Dki′中的分水嶺所在的位置；由這些分水嶺即把圖像Dki′分割成若干個子區(qū)域；步驟(4.4).對分水嶺算法得到的圖像Dki′中的每個子區(qū)域，分別計算其與下方相鄰層面器官分割區(qū)域Oki-1的面積重合比率 把重合比率低于50％的區(qū)域視為錯誤分割區(qū)域，重新標記為背景區(qū)域，從而將其在距離圖像Dki′中除去。
全文摘要
本發(fā)明屬于醫(yī)學圖像中腹部器官分割技術領域。其特征在于，該方法依次含有以下步驟采用基于歸一化互信息的配準方法將標準圖譜向通過CT或核磁掃描得到的個體圖像作整體配準；用相似性測度對每個感興趣器官進行單獨配準；模糊連接分割和器官形狀修正；其中，模糊連接分割又分為先算出曲線擬和后的灰度直方圖曲線，并視其為目標圖像中灰度概率密度的函數(shù)；再把具有最大區(qū)域灰度概率密度的點作為種子點，將直方圖擬和曲線以及種子作為初始參數(shù)，用模糊連接算法得到模糊連接圖像，再用設定的閾值把目標像素和背景像素分開。本發(fā)明具有可同時進行目標圖像中腹部多器官自動分割，實施簡單的優(yōu)點。
文檔編號A61B5/055GK101013503SQ20071006306
公開日2007年8月8日 申請日期2007年1月26日 優(yōu)先權日2007年1月26日
發(fā)明者白凈, 周永新, 王洪凱, 劉加成, 張永紅, 張菊鵬 申請人:清華大學