專利名稱：一種對管狀結(jié)構(gòu)內(nèi)層分割的方法
技術(shù)領(lǐng)域：
本發(fā)明涉及計算機(jī)圖像處理領(lǐng)域，特別涉及一種基于模糊聚類和活 動輪廓的結(jié)合縱向和橫向截平面的管狀結(jié)構(gòu)內(nèi)層分割的方法。
背景技術(shù)：
圖像分割是指將圖像中若干感興趣的區(qū)域提取出來，感興趣的區(qū)域 通常稱為前景，其它區(qū)域稱為背景。在分割的基礎(chǔ)上才可能對目標(biāo)進(jìn)行 定量分析或者識別，從而使更高層的分析和理解成為可能。
九十年代以來，隨著技術(shù)的發(fā)展，可以獲得更高分辨率的圖像，基 于形變模型的方法也開始大量應(yīng)用。形變模型是一種尋找圖像中物體邊
界的方法。給定初始輪廓，然后通過圖像信息及先驗知識定義外力，推 動輪廓到達(dá)物體的邊界。形變模型分為兩大類參數(shù)形變模型和幾何形 變模型。
參數(shù)化形變模型最早可以追溯到Kass等人提出的活動輪廓模型。分 割過程就是模型在外力和內(nèi)力的作用下向物體邊緣靠近。內(nèi)力控制曲線 的形狀，外力推動曲線運動。然而傳統(tǒng)的活動輪廓模型對初始輪廓的敏 感性及有限的捕捉范圍c對此，許多學(xué)者提出了多種外力的定義。M.Kass
和D. Terzopoulos使用高斯函數(shù)對圖像進(jìn)行平滑，雖然擴(kuò)大了捕捉范圍， 但是隨著C7的增加，可能會使邊緣平滑的太厲害，從而找不到真實的邊 緣。L.D.Cohen等提出了 balloon模型，在外力中增加了壓力(pressure force)項，該模型改善了活動輪廓對初始輪廓的敏感性，并且能夠跨越 偽邊緣點，但是需要控制壓力的大小及方向。而且該模型無法收斂到深 度凹陷部分。L.D.Cohen和I.Cohen提出了距離力可以擴(kuò)大捕捉范圍，但 是距離力使模型上的點靠近最近的邊界點，這使得輪廓無法收斂到深度凹陷處。
1998年Xu等提出了梯度矢量流來解決不能收斂到深度凹陷處的問 題，他運用向量擴(kuò)散公式將梯度擴(kuò)散到遠(yuǎn)離邊緣的區(qū)域中，于是產(chǎn)生了
另外一個力場，被稱為梯度矢量流(gradient vector flow, GVF)場。GVF 改善了邊緣凹陷問題。雖然GVF解決了活動輪廓捕捉范圍小的問題，但 是如果初始輪廓位置不正確的話，仍然有可能導(dǎo)致活動輪廓收斂到錯位 的位置。這導(dǎo)致了基于活動輪廓的方法都需要人工設(shè)定初始輪廓，然而 如果圖像的數(shù)據(jù)量很大，人工設(shè)定初始輪廓即費力又費力，重復(fù)性又差， 因此必須使用自動化的方法進(jìn)行分割。

發(fā)明內(nèi)容
本發(fā)明的目的是提供一種自動的圖像分割方法，能夠?qū)軤罱Y(jié)構(gòu)內(nèi) 層進(jìn)行準(zhǔn)確、自動分割的方法。
為實現(xiàn)上述目的，本發(fā)明的一種對管狀結(jié)構(gòu)內(nèi)層分割的方法，結(jié)合 縱向和橫向截平面的管狀結(jié)構(gòu)內(nèi)層自動分割的步驟包括
步驟a:對管狀結(jié)構(gòu)三維序列數(shù)據(jù)進(jìn)行各向異性擴(kuò)散濾波，用于除 去隨機(jī)噪聲；
步驟b:從管狀結(jié)構(gòu)三維序列數(shù)據(jù)中提取縱向截平面，并在縱向截 平面上利用模糊聚類提取管狀結(jié)構(gòu)內(nèi)層的輪廓線；
步驟C:把縱向截平面提取的輪廓線變換到橫向截平面，每一個橫 向截平面得到對應(yīng)的輪廓點；
步驟d:對橫向截平面上的輪廓點進(jìn)行樣條插值，得到一條封閉的
樣條輪廓；
步驟e:利用小波變換產(chǎn)生橫向截平面的多尺度圖像，用來在多尺 度圖像中使用活動輪廓，多尺度圖像包含三個尺度，即第一尺度、第二 尺度和第三尺度；
步驟f:把樣條輪廓變換到不同尺度，用來作為多尺度圖像中活動
輪廓的初始輪廓；
步驟g:計算橫向截平面的多尺度圖像的梯度矢量流場；
步驟h:在橫向截平面的多尺度圖像中利用活動輪廓進(jìn)化輪廓，活動輪廓的外力使用梯度矢量流場；
步驟h將多尺度圖像中活動輪廓進(jìn)化得到的輪廓通過變換及加權(quán)， 得到管狀結(jié)構(gòu)的最終內(nèi)層輪廓。
優(yōu)選的實施例，在縱向截平面提取內(nèi)層輪廓線的步驟包括
步驟bl:從管狀結(jié)構(gòu)三維序列數(shù)據(jù)中等角度的方向上提取 乂* [3,6]個縱向截平面，W為需要處理的縱向截平面的個數(shù)。
步驟b2:移除縱向截平面中不包含有用信息的導(dǎo)管區(qū)域，得到兩個 縱向子截平面；
步驟b3:使用"*"， " = {3,5,7,9}的滑動窗口，滑動窗口中心點遍歷縱 向子截平面的每一個像素點，并統(tǒng)計窗口內(nèi)像素灰度的均值和方差，并 以統(tǒng)計得到的均值和方差作為特征進(jìn)行模糊聚類，"為滑動窗口的邊長;
步驟b4:利用模糊聚類把每一個縱向子截平面分為兩類， 一類為內(nèi) 層區(qū)域，另一類為管壁及其它區(qū)域，從這兩類的邊界可提取管狀結(jié)構(gòu)在 縱向截平面的內(nèi)層輪廓線。
優(yōu)選的實施例，其中，提取最終內(nèi)層輪廓的步驟包括
步驟il:將不同尺度的活動輪廓進(jìn)化得到的輪廓變換到第一尺度(即 原始圖像所在的尺度)，得到不同尺度的內(nèi)層輪廓；
步驟i2:對不同尺度的內(nèi)層輪廓分別進(jìn)行形態(tài)學(xué)操作，填充輪廓， 得到內(nèi)層區(qū)域；
步驟i3:對不同尺度的內(nèi)層區(qū)域進(jìn)行距離變換，得到內(nèi)層區(qū)域距離
步驟i4:對內(nèi)層區(qū)域距離圖A'^"3進(jìn)行加權(quán)，通過以下公式進(jìn)行 加權(quán)
其中M、 ^、 ^為加權(quán)因子，加權(quán)因子的取值范圍為 0《A <1'0 < 2《i'o《巧《1， < i十^ =1 ，并設(shè)定 一 個閨值 T， ^(幼1'効2'%}0<1，最終內(nèi)層區(qū)域可以通過公式D加f =">7>尋到，
大于該闊值的區(qū)域即為管狀結(jié)構(gòu)的最終內(nèi)層區(qū)域，小于該閾值的區(qū)域被 舍棄，^>"/為最終內(nèi)層區(qū)域；最終內(nèi)層輪廓即為最終內(nèi)層區(qū)域的邊界。
本發(fā)明的有益效果本發(fā)明利用模糊聚類在縱向截平面上提取管狀結(jié)構(gòu)內(nèi)層的輪廓線，然后把縱向截平面的輪廓線變換到橫向截平面，進(jìn) 行輪廓的初始化。不僅避免了手工初始化管狀結(jié)構(gòu)內(nèi)層的輪廓，而且大 大提高了初始化管狀結(jié)構(gòu)內(nèi)層輪廓的速度?；顒虞喞惴ㄒM可能減少 噪聲的干擾，然而圖像噪聲的存在，導(dǎo)致活動輪廓不能得到精確的內(nèi)層 輪廓。為了減少噪聲對活動輪廓的影響，首先，使用各向異性濾波去除 隨機(jī)噪聲，各向異性擴(kuò)散濾波能在去除隨機(jī)噪聲的同時還能保持邊緣信 息并能增強邊緣；其次，利用小波變換產(chǎn)生多尺度圖像，并在不同尺度 上利用活動輪廓進(jìn)化輪廓。最后，把不同尺度上得到的輪廓加權(quán)變換到 原始圖像尺度，得到管狀結(jié)構(gòu)內(nèi)層的最終輪廓。本發(fā)明利用了模糊聚類 和活動輪廓在多尺度上提取管狀結(jié)構(gòu)的內(nèi)層輪廓，在保持了抗噪性能好 的情況下，實現(xiàn)了初始化的自動化，從而實現(xiàn)了整個分割方法的自動化。


此處所說明的附圖用來提供對本發(fā)明的進(jìn)一步理解，構(gòu)成本申請的 一部分，并不構(gòu)成對本發(fā)明的不當(dāng)限定，在附圖中
圖1是本發(fā)明的管狀結(jié)構(gòu)內(nèi)層分割方法的實現(xiàn)流程圖2是本發(fā)明實施例中提取縱向截平面的示意圖3a—圖3f是本發(fā)明實施例中在縱向截平面上提取管狀結(jié)構(gòu)內(nèi)層
輪廓線的示意圖4是本發(fā)明實施例中縱向截平面管狀結(jié)構(gòu)內(nèi)層輪廓點與橫向截平
面對應(yīng)關(guān)系的示意圖5a—圖5c是本發(fā)明實施例中樣條插值的示意圖6a—圖6c是本發(fā)明實施例中多尺度圖像中活動輪廓進(jìn)化示意圖7是本發(fā)明實施例中提取管狀結(jié)構(gòu)最終內(nèi)層輪廓示意圖。
具體實施例方式
下面結(jié)合附圖詳細(xì)說明本發(fā)明技術(shù)方案中所涉及的各個細(xì)節(jié)問題。 所描述的實施例僅旨在便于對本發(fā)明的理解，而對其不起任何限定作用。
作為一種具體的實現(xiàn)方案，使用€++語言實現(xiàn)了本發(fā)明所描述的算 法，所有的實現(xiàn)都是在一臺酷睿2雙核2.2G， 2G內(nèi)存，操作系統(tǒng)為Windows XP的計算機(jī)上完成的。圖1示出了可用于實施在此公開的方
法的流程圖。該方法包括
步驟100:讀入三維序列數(shù)據(jù)；
步驟101:各向異性擴(kuò)散濾波如下所述圖像中常含有隨機(jī)噪聲， 因此圖像預(yù)處理是圖像分析與識別必要的前期工作。為了減少噪聲對分 割算法的影響，采用各向異性擴(kuò)散濾波來去除管狀結(jié)構(gòu)三維序列數(shù)據(jù)的 隨機(jī)噪聲，同時還能保持邊緣信息和增強邊緣位置。各向異性擴(kuò)散濾波
(如下所述P. Perona and J. Malik, "Scale-space and edge detection using anisotropic diffusion," IEEE Trans. Pattern Anal. Machine Intell., vol, 12, no. 7 pp. 629~639， July 1990)最早是由Perona和Malik提出的，可以用下面 的微分模型表示
其中x是圖像的橫坐標(biāo)，y是圖像的縱坐標(biāo)，/是時間變量，力'v是散 度算子，v/是圖像的辨度。擴(kuò)散系數(shù)是關(guān)于梯度幅值的非負(fù)單調(diào)遞減函
數(shù)g(夂滿足-)=1,^*)=()，這樣可以保證區(qū)域內(nèi)部(梯度小)，擴(kuò)
散系數(shù)大；在邊緣(梯度大)，擴(kuò)散系數(shù)小，從而起到邊緣增強的作用。 由于Perona和Malik提出的擴(kuò)散系數(shù)比較大，可能會消除一些不明顯而 又重要的邊緣細(xì)節(jié)，本發(fā)明采用基于Tukey的Biweight Estimator誤差模 型(如下所述Black M. J.， Snpiro G.. Marimont D. H., Heeger D" "Robust anisotrupic diffbsiun". IEEE Trims. On Image Pmc., Vol. 7.No 3， pp. 421-432, March 1998)的邊緣終止函數(shù)為
g(V/,CT)=
1^^j作為擴(kuò)散系數(shù)，其中"是梯度閾值，
其它
基于Tukey函數(shù)比Perona-Malik函數(shù)遞減程度更快，當(dāng)梯度超過cr時， 擴(kuò)散系數(shù)為0，即模型不再擴(kuò)散，從而避免了邊緣弱化問題。
步驟102:從管狀結(jié)構(gòu)三維序列數(shù)據(jù)中提取縱向截平面，并在縱向 截平面上利用模糊聚類提取管狀結(jié)構(gòu)內(nèi)層的輪廓線如下所述模糊C均 值聚類(FCM)算法是Bezkek(如下所述J.C. Bezdek,Pattern recognition with fuzzy objective fUnction algorithms,Plenum Press, New York; 1981)于
81981年提出的。其基本思想是".,/ = 1,2,...，^}是^個樣本組成的樣本集 合I， C為預(yù)定的類別數(shù)目，C/,/ = l,2,...，C為每個聚類的中心， 是第/ 個樣本對第y個樣本的隸屬度函數(shù)。目標(biāo)準(zhǔn)則函數(shù)可以寫為
人-SZ《k-。1 ，其中/ >1是一個可以控制聚類結(jié)果的模糊程度
的常數(shù)。模糊c均值算法是通過對目標(biāo)準(zhǔn)則函數(shù)的迭代來獲取對數(shù)據(jù)集
的模糊分類，即迭代
_。
2>;
使目標(biāo)準(zhǔn)則函數(shù)收斂到一個局部
'水
極小值或鞍點，得到x的一個最優(yōu)模糊c劃分。 是第/個樣本對第.
樣本的隸屬度函數(shù)，C為聚類數(shù)目，。為第/個聚類中心，q為第A個聚 類中心，x,為第/個樣本，附為控制聚類結(jié)構(gòu)的模糊程度的常數(shù)，且附>1， W為樣本的數(shù)目。
由于噪聲的影響，僅利用像素的灰度信息對圖像進(jìn)行聚類，聚類結(jié)果 對噪聲較為敏感。為此，使用一個fl、 " = {3,5，7,9}滑動窗口，其中"為滑 動窗口的邊長，滑動窗口的中心點遍歷縱向截平面的每一個像素，并統(tǒng)計 窗口內(nèi)的均值和方差作為模糊聚類的特征A，即《={5, }， ^為窗口內(nèi) 像素的均值，<^為窗口內(nèi)像素的方差。以滑動窗口內(nèi)的均值和方差作為特 征進(jìn)行模糊聚類，不僅利用了像素的灰度信息，還利用了窗口內(nèi)像素的方 差信息，對噪聲有一定的抑制能力。
活動輪廓分割圖像需要初始輪廓。然而管狀結(jié)構(gòu)三維序列數(shù)據(jù)有上百 張乃至上千張的圖像，如果每張都手動描繪初始輪廓，這樣工作量大，耗 時多，而且可重復(fù)性差。有些半自動的方法需要初始化第一張圖像的輪廓， 然后利用上一張圖像的分割結(jié)果作為下一張圖像的初始輪廓，如果相鄰兩 張圖像的輪廓變化較大，這樣會導(dǎo)致初始化位置不正確，可能出現(xiàn)分割錯 誤。為此，本發(fā)明利用模糊聚類在縱向截平面上提取管狀結(jié)構(gòu)內(nèi)層的輪廓 線，然后變換到橫向截平面初始輪廓。每一張圖像的初始輪廓都是獨立的， 每一張圖像的初始輪廓不會受相鄰圖像初始輪廓的影響，包括以下步驟步驟102a:請參考圖2提取縱佝截平面的示意圖，提取三維序列數(shù) 據(jù)的縱向截平面是從管狀結(jié)構(gòu)三維序列數(shù)據(jù)20中等角度的方向上提取 iV，臉f3，6]個縱向截平面21,22， iV為需要處理的縱向截平面的個數(shù)。
步驟102b:移除縱向截平面的導(dǎo)管區(qū)域如下所述由于縱向截平面 21， 22的圖像中包括三個區(qū)域?qū)Ч軈^(qū)域，內(nèi)層區(qū)域、管壁和其它區(qū)域。 導(dǎo)管區(qū)域不包含任何有用的信息，可以移除導(dǎo)管區(qū)域?？v向截平面22 移除導(dǎo)管區(qū)域后得到左、右兩個縱向子截平面30。
步驟102c:請參考圖3a—圖3f本發(fā)明實施例中在縱向截平面22上 提取管狀結(jié)構(gòu)內(nèi)層輪廓線的示意圖；其中如圖3a示出 個縱向截平面， 在縱向截平面上進(jìn)行模糊聚類如下所述使用^a， 。 = {3，5，7，9}的滑動窗 口，其中a為滑動窗口的邊長；如圖3b示出移除導(dǎo)管區(qū)域后得到的兩個 縱向子截平面，滑動窗口的中心點遍歷縱向子截平面30的每一個像素， 并統(tǒng)計窗口內(nèi)像素的均值和方差。并以統(tǒng)計得到的均值和方差作為特征 進(jìn)行模糊聚類。每個縱向子截平面聚為兩類， 一類為內(nèi)層區(qū)域， 一類為 管壁和其它區(qū)域。如圖3c示出了模糊聚類的結(jié)果。由于斑點噪聲的影響， 雖然使用了均值和方差作為特征進(jìn)行聚類，但是聚類結(jié)果中可能還會出 現(xiàn)噪聲點。因此使用形態(tài)學(xué)操作消除噪聲點和不規(guī)則的邊緣，聚類后的 圖像31經(jīng)過形態(tài)學(xué)操作得到縱向截平面管狀結(jié)構(gòu)內(nèi)層的分割結(jié)果32， 如圖3d示出了形態(tài)學(xué)操作后的圖像。
步驟102d:提取內(nèi)層輪廓線如下所述利用模糊聚類把每一個縱向 子截平面分為兩類， 一類為內(nèi)層區(qū)域，另一類為管壁及其它區(qū)域，從這 兩類的邊界可提取管狀結(jié)構(gòu)在縱向截平面的內(nèi)層輪廓線。續(xù)請參考圖 3d，對于左子截平面321，從右到左掃描圖像，當(dāng)像素從黑色變?yōu)榘咨?時，該白色點即為輪廓點，然后跳轉(zhuǎn)到下一行掃描。對于右子截平面322， 從左到右進(jìn)行掃描。兩個子截平面都處理完畢后，得到縱向截平面上的 內(nèi)層輪廓線33，如圖3e示出提取的管狀結(jié)構(gòu)內(nèi)層輪廓線。輪廓線疊加 到原始縱向截平面圖像34，如圖3f示出管狀結(jié)構(gòu)內(nèi)層輪廓線疊加到原 始縱向截平面。
步驟103:把縱向截平面提取的輪廓線變換到橫向截平面，每一個 橫向截平面得到對應(yīng)的輪廓點；如圖4示出本發(fā)明實施例中縱向截平面管狀結(jié)構(gòu)內(nèi)層輪廓點和橫向截平面對應(yīng)關(guān)系的示意圖；將縱向截平面內(nèi) 層輪廓線變換到橫向截平面如下所述為了能在橫向截平面使用活動輪 廓，需要在橫向截平面上初始輪廓，因此，需要把縱向截平面提取的輪
廓線變換到橫向截平面?？v向截平面提取的輪廓線40與橫向截平面42
有唯一的對應(yīng)關(guān)系，兩個"☆"標(biāo)志為同一輪廓點在縱向截平面和橫向
截平面中的位置41。本發(fā)明提取了管狀結(jié)構(gòu)三維序列數(shù)據(jù)中等角度的 N-5個縱向截平面，每一個縱向截平面能提取兩條內(nèi)層輪廓線。把五個 縱向截平面的輪廓線變換到橫向截平面，每一個橫向截平面會得到對應(yīng) 的十個輪廓點50。
步驟104:如圖5a—圖5c示出本發(fā)明實施例中樣條插值的示意圖； 其中圖5a示出橫向截平面的輪廓點，對橫向截平面上的輪廓點進(jìn)行樣 條插值，得到一條封閉的樣條輪廓；如下所述；縱向截平面的輪廓線變 換到橫向截平面后，只有十個輪廓點50，初始輪廓需要的是一條封閉的 曲線，而不是獨立的點。因此，需要對輪廓點進(jìn)行樣條插值。三次樣條 在靈活性和計算速度之間提供了一個合理的折衷方案，與低次多項式相 比，三次樣條在模擬任意曲線時顯得更靈活，與更高次多項式相比，三 次樣條只需較少的計算和存儲空間。因此，初始輪廓使用Cardinal三次
樣條(如下所述Donald Heam ， M Pauline Baker . Computer Graphics[M]. Prentice Hall Press Company. 1998)進(jìn)行插值得到的樣條輪 廓。Cardinal樣條完全由連續(xù)4個控制點給出，中間2個點是曲線的端 點，第1個點和第4個點用來計算斜率。插值公式可以由以下多項式表 示
尸(")=A-iO"3 +2ra2 -卯)+凡[(2-*3 +(s-3)w2 +1]
其中^(l-0/2，,為張力參數(shù)，"為點yt和點^ + l間的相對偏移量，
， ，)為相對偏移量是"時，插值點的位置。A,i^為曲線的端點，尸w，尸w 為連續(xù)四個點的第1個點和第4個點。輪廓點使用上述公式進(jìn)行插值得 到樣條輪廓51，如圖5b示出Cardinal三次樣條插值得到的樣條輪廓。 把樣條輪廓疊加到原始橫向截平面圖像52，如圖5c示出樣條插值的輪 廓疊加到原始橫向截平面。
ii步驟105:利用小波變換產(chǎn)生橫向截平面的多尺度圖像，用來在多 尺度圖像中使用活動輪廓；小波變換得到多尺度圖像如下所述輪廓進(jìn) 化要盡量減少噪聲的干擾，并同時盡量多保留感興趣的邊緣。但是兩者 往往不能同時兼顧。多尺度分析提供了有效的解決方法。小尺度時圖像 比較精細(xì)，但是噪聲也比較嚴(yán)重，大尺度時圖像比較粗略，但噪聲得到 抑制。因此，可以使用多尺度的方法提取管狀結(jié)構(gòu)的內(nèi)層輪廓。使用小 波變換生成^層小波系數(shù)^X，2。將原始橫向截平面圖像記為第一尺
度圖像42， ^為第二尺度圖像61， 4為第三尺度圖像63，則第一尺度 圖像42、第二尺度圖像61、第三尺度圖像63可以看成是由精細(xì)到粗略 的多尺度圖像。
步驟106:把樣條輪廓變換到不同尺度，用來作為多尺度圖像中活 動輪廓的初始輪廓；將樣條輪廓變換到不同尺度如下所述步驟104得 到了一個樣條輪廓51，這個樣條輪廓51只能作為第一尺度圖像42的初 始輪廓，第二尺度圖像6K第三尺度圖像63的初始輪廓還沒有得到。 因此，需要把樣條輪廓51變換到第二尺度和第三尺度。將樣條輪廓51 變換到第二尺度，得到第二尺度圖像的初始輪廓60，將樣條輪廓51變 換到第三尺度，得到第三尺度圖像的初始輪廓62。
步驟107:計算橫向截平面的多尺度圖像的梯度矢量流場。梯度矢 量流場(GVF Field)(如下所述XU C Y, P rince J L. Snakes, shapes, and
gradient vector flow [ J ]. IEEE T ransactions on Imaging Processing, 1995)
可以表示為矢量場v(;c，力-["(x，力,v(:c,力]，并且滿足以下能量最小化函數(shù):
/是邊緣圖像，V/是邊緣圖像的梯度，v為矢量場，",v分別為矢量場的 第一和第二個分量，^, ,^,^是^分別對x,少求一階偏導(dǎo)，s為能量最 小化函數(shù)，/v為平衡能量最小化函數(shù)中前后兩項權(quán)重的系數(shù)，A隨著噪 聲的增加而增大。依據(jù)變分法，梯度矢量流場可以通過求解下列Euler 方程得到-
其中，/;,/y為邊緣圖像/對^的偏導(dǎo)，A是拉普拉斯算子。以梯度矢量
盧-(V-/》|y/f=o流場作為外力極大地增加了活動輪廓的捕獲范圍。每一個尺度的圖像都 需要計算對應(yīng)尺度的梯度向量流場。
步驟108:在橫向截平面的多尺度圖像中利用活動輪廓進(jìn)化輪廓；
如圖6a—6c示出了多尺度圖像中活動輪廓進(jìn)化示意圖，在多尺度圖像中
利用活動輪廓進(jìn)化輪廓如下所述
活動輪廓是由Kass等人提出的(如下所述M.Kass,A. Witkin,and D. Terzopoulos, "Snakes: Active contourmodels,，, Int. J. Comput. Vis., vol. 1， pp. 321-331, 1987 )，活動輪廓定義為一條參數(shù)化曲線 X(" = (;^),3^)),K
，曲線在圖像中移動并滿足下面最小化能量函數(shù)
其中，￡為最小化能量函數(shù)，義("為參數(shù)化曲線，"為加權(quán)因子，々為 加權(quán)因子，￡ew(X(a-))為與圖像有關(guān)的外力，X'(s)為關(guān)于參數(shù)s的一階導(dǎo) 數(shù)，JT(s)為關(guān)于參數(shù)s的二階導(dǎo)數(shù)，最小化能量函數(shù)￡必須滿足以下Euler 方程
其中;r(s)為關(guān)于參數(shù)s的二階導(dǎo)數(shù)，r"(s)為關(guān)于參數(shù)s的四階導(dǎo)數(shù)，
V為梯度算子，取V&w-v(;c,力，v(;c,力為步驟107計算的梯度矢量流場。 參數(shù)化曲線的初始位置即為活動輪廓的初始輪廓，步驟106可以得 到不同尺度活動輪廓的初始輪廓。通過求解上述的Euler方程，可以得 到輪廓進(jìn)化的結(jié)果。
由于噪聲的干擾，單純的應(yīng)用活動輪廓并不能得到精確的輪廓。因 此，本發(fā)明在三個不同尺度上使用活動輪廓進(jìn)化輪廓。三個尺度的輪廓 進(jìn)化可以同時進(jìn)行，這樣可以使用并行程序?qū)崿F(xiàn)活動輪廓進(jìn)化，減少了 輪廓進(jìn)化時間。第一尺度圖像42和第一尺度的初始輪廓51使用活動輪 廓進(jìn)化得到的輪廓611，如圖6a示出第一尺度活動輪廓進(jìn)化的示意圖、 第二尺度圖像61和第二尺度的初始輪廓60使用活動輪廓進(jìn)化得到的輪 廓621，如圖6b示出第二尺度活動輪廓進(jìn)化的示意圖、第三尺度圖像63 和第三尺度的初始輪廓62使用活動輪廓進(jìn)化得到的輪廓631，如圖6c 示出第三尺度活動輪廓進(jìn)化的示意圖。步驟109:如圖7示出了提取管狀結(jié)構(gòu)最終內(nèi)層輪廓示意圖，將多 尺度圖像中活動輪廓進(jìn)化得到的輪廓通過變換及加權(quán)，得到管狀結(jié)構(gòu)的 最終內(nèi)層輪廓如下所述不同尺度的圖像有不同的特征，小尺度圖像比 較精細(xì)，大尺度圖像比較粗略。因此，本發(fā)明結(jié)合不同尺度的分割結(jié)果， 得到最終內(nèi)層的輪廓。首先，將不同尺度的活動輪廓進(jìn)化得到的輪廓變 換到第一尺度(即原始圖像所在的尺度)，由于不同尺度圖像的大小不一 樣，為了結(jié)合不同尺度活動輪廓進(jìn)化得到的輪廓，需要把第二尺度和第 三尺度活動輪廓進(jìn)化得到的輪廓變換到第一尺度上，第一尺度活動輪廓 進(jìn)化得到的輪廓不需要變換，因此第一尺度內(nèi)層輪廓611即為第一尺度
活動輪廓進(jìn)化得到的輪廓611，第二尺度活動輪廓進(jìn)化得到的輪廓621 變換到第一尺度，得到第二尺度內(nèi)層輪廓622;第三尺度活動輪廓進(jìn)化 得到的輪廓631變換到第一尺度，得到第三尺度內(nèi)層輪廓632。
其次，對不同尺度的內(nèi)層輪廓分別進(jìn)行形態(tài)學(xué)操作，填充輪廓，得 到內(nèi)層區(qū)域。第一尺度內(nèi)層輪廓6U進(jìn)行區(qū)域填充，得到第一尺度內(nèi)層 區(qū)域72;第二尺度內(nèi)層輪廓622進(jìn)行區(qū)域填充，得到第二尺度內(nèi)層區(qū)域 71;第三尺度內(nèi)層輪廓632進(jìn)行區(qū)域填充，得到第三尺度內(nèi)層區(qū)域70;
然后，對不同尺度的內(nèi)層區(qū)域進(jìn)行距離變換，得到內(nèi)層區(qū)域距離圖 A,,' = l,2,3， D,為第一尺度內(nèi)層區(qū)域距離圖75， D2為第二尺度內(nèi)層區(qū)域 距離圖74， ￡>3為第三尺度內(nèi)層區(qū)域距離圖73。
最后，對內(nèi)層區(qū)域距離圖A,^1,2，3進(jìn)行加權(quán)，通過以下公式進(jìn)行加
權(quán)
<formula>formula see original document page 14</formula> 其中w, 、 《2 、約為加權(quán)因子，加權(quán)因子的取值范圍為 0 ,《，0 一《，0 一《，+ 2 + w3 = 1 ，并設(shè)定 一 個閾值 r, min^，"2,約)^r^1，最終內(nèi)層區(qū)域可以通過公式1%"/=1>>7得到，
大于該閾值的區(qū)域即為管狀結(jié)構(gòu)的最終內(nèi)層區(qū)域，小于該閾值的區(qū)域被 舍棄，/為最終內(nèi)層區(qū)域。通過閾值判斷，得到管狀結(jié)構(gòu)最終內(nèi)層區(qū)
域77。最終管狀結(jié)構(gòu)的內(nèi)層輪廓78即為管狀結(jié)構(gòu)最終內(nèi)層區(qū)域77的邊 界。
本發(fā)明采用的是從不同縱向截平面提取內(nèi)層輪廓線，并把縱向截平面的輪廓線變換到橫向截平面來初始化活動輪廓的輪廓，而不是使用上 一張圖像的輪廓作為下一張圖像的初始輪廓，這樣避免了初始化錯誤， 提高了分割的精度。由于每張圖像都有一個初始輪廓，可以使用并行程 序同時對多張圖像的輪廓進(jìn)行進(jìn)化。本發(fā)明還采用了多尺度分割管狀結(jié) 構(gòu)內(nèi)層來提高分割精度。
以上所述，僅為本發(fā)明中的具體實施方式
，但本發(fā)明的保護(hù)范圍并 不局限于此，任何熟悉該技術(shù)的人在本發(fā)明所揭露的技術(shù)范圍內(nèi)，可理 解想到的變換或替換，都應(yīng)涵蓋在本發(fā)明的包含范圍之內(nèi)，因此，本發(fā) 明的保護(hù)范圍應(yīng)該以權(quán)利要求書的保護(hù)范圍為準(zhǔn)。
權(quán)利要求
1.一種對管狀結(jié)構(gòu)內(nèi)層分割的方法，其特征在于，包括步驟步驟a對管狀結(jié)構(gòu)三維序列數(shù)據(jù)進(jìn)行各向異性擴(kuò)散濾波，用于除去隨機(jī)噪聲；步驟b從管狀結(jié)構(gòu)三維序列數(shù)據(jù)中提取縱向截平面，并在縱向截平面上利用模糊聚類提取管狀結(jié)構(gòu)內(nèi)層的輪廓線；步驟c把縱向截平面提取的輪廓線變換到橫向截平面，每一個橫向截平面得到對應(yīng)的輪廓點；步驟d對橫向截平面上的輪廓點進(jìn)行樣條插值，得到一條封閉的樣條輪廓；步驟e利用小波變換產(chǎn)生橫向截平面的多尺度圖像，用來在多尺度圖像中使用活動輪廓，多尺度圖像包含第一尺度、第二尺度和第三尺度；步驟f把樣條輪廓變換到不同尺度，用來作為多尺度圖像中活動輪廓的初始輪廓；步驟g計算橫向截平面的多尺度圖像的梯度矢量流場；步驟h在橫向截平面的多尺度圖像中利用活動輪廓進(jìn)化輪廓，活動輪廓的外力使用梯度矢量流場；步驟i將多尺度圖像中活動輪廓進(jìn)化得到的輪廓通過變換及加權(quán)，得到管狀結(jié)構(gòu)的最終內(nèi)層輪廓。
2. 按照權(quán)利要求l所述方法，其特征在于，在縱向截平面提取內(nèi)層輪廓線的步驟包括-步驟bl:從管狀結(jié)構(gòu)三維序列數(shù)據(jù)中等角度的方向上提取 iV,A^ [3,6]個縱向截平面，W為需要處理的縱向截平面的個數(shù)。步驟b2:移除縱向截平面中不包含有用信息的導(dǎo)管區(qū)域，得到兩個縱向子截平面；步驟b3:使用Wfl， " = {3,5,7,9}的滑動窗口，滑動窗口中心點遍歷縱 向子截平面的每一個像素點，并統(tǒng)計窗口內(nèi)像素灰度的均值和方差，并 以統(tǒng)計得到的均值和方差作為特征進(jìn)行模糊聚類，"為滑動窗口的邊長;步驟b4:利用模糊聚類把每一個縱向子截平面分為兩類， 一類為內(nèi) 層區(qū)域，另一類為管壁及其它區(qū)域，從這兩類的邊界可提取管狀結(jié)構(gòu)在 縱向截平面的內(nèi)層輪廓線。
3.按照權(quán)利要求1所述方法，其特征在于，所述提取最終內(nèi)層輪廓 的步驟包括步驟il:將不同尺度的活動輪廓進(jìn)化得到的輪廓變換到第一尺度(即 原始圖像所在的尺度)，得到不同尺度的內(nèi)層輪廓；步驟i2:對不同尺度的內(nèi)層輪廓分別進(jìn)行形態(tài)學(xué)操作，填充輪廓， 得到內(nèi)層區(qū)域；步驟i3:對不同尺度的內(nèi)層區(qū)域進(jìn)行距離變換，得到內(nèi)層區(qū)域距離 圖q.'i =1,2,3:步驟i4:對內(nèi)層區(qū)域距離圖A'i^'2'3進(jìn)行加權(quán)，通過以下公式進(jìn)行 加權(quán)其中A 、 《2 、約為加權(quán)因子，加權(quán)因子的取值范圍為 0 <巧<1, 0 <<2>2 <1, 0 <巧<1，i+ 2+iJ>3=1 ， 并設(shè)定 一 個閾值r，腿(巧,巧'巧)〈7^1，最終內(nèi)層區(qū)域可以通過公式/)—廣i)〉:r得到，大于該閾值的區(qū)域即為管狀結(jié)構(gòu)的最終內(nèi)層區(qū)域，小于該閾值的區(qū)域被 舍棄，z)細(xì),為最終內(nèi)層區(qū)域；最終內(nèi)層輪廓即為最終內(nèi)層區(qū)域的邊界。
全文摘要
本發(fā)明公開了一種對管狀結(jié)構(gòu)內(nèi)層分割的方法，該方法首先提取管狀結(jié)構(gòu)三維序列數(shù)據(jù)的縱向截平面，并在縱向截平面提取管狀結(jié)構(gòu)的內(nèi)層輪廓線，然后把縱向截平面的輪廓線變換到橫向截平面，接著把橫向截平面的輪廓點樣條插值得到樣條輪廓，利用小波變換產(chǎn)生橫向截平面的多尺度圖像，并把樣條輪廓變換到不同尺度作為活動輪廓的初始輪廓，然后進(jìn)化輪廓，最后結(jié)合不同尺度活動輪廓的分割結(jié)果，得到管狀結(jié)構(gòu)的最終內(nèi)層輪廓。應(yīng)用本發(fā)明的技術(shù)方案減少了初始輪廓的時間。為提高分割精度，該方法在多尺度上使用活動輪廓進(jìn)化輪廓。
文檔編號G06T7/00GK101599176SQ20091008735
公開日2009年12月9日 申請日期2009年6月17日 優(yōu)先權(quán)日2009年6月17日
發(fā)明者何暉光, 易建華 申請人:中國科學(xué)院自動化研究所