萬小萍 顧 勇（重慶大學(xué)生物工程學(xué)院，重慶 400044）

2（重慶警官職業(yè)學(xué)院，重慶 400039）

引言

圖像分割是圖像處理過程中具有挑戰(zhàn)性的關(guān)鍵步驟，其目的在于將一幅圖像按照一定的特征（如顏色、紋理等）分成若干有意義的互不交疊的同質(zhì)區(qū)。在各種分割方法中，由Kass等在1988年提出的可變形模型［1］，由于其能很好的分割出興趣目標(biāo)的形狀，近年來成為圖像分割中最常用的方法之一?？勺冃文Ｐ驮谠S多場合都能得到滿意的結(jié)果，但是在有些情況下，如圖像中存在大量的噪聲、目標(biāo)物體有缺損或含有錯誤信息時，大部分不含先驗知識的算法都很難得到有意義的分割結(jié)果。

此時要想得到滿意的分割結(jié)果，其中一個比較好的解決辦法是在分割過程中融入形狀先驗知識。目前，已有大量相關(guān)文獻報道，其中基于水平集［2］，Mumford-Shah（M-S）模型［3］和先驗形狀［4］等具有一定的代表性。2002年，Cremers等介紹了一種將形狀先驗和M-S模型結(jié)合的統(tǒng)計模型［5］，這種模型將曲線顯性地表示為一條封閉的樣條曲線。2003年，Tsai等提出的形狀分割模型［6］，基于Leventon的形狀模型［7］和Chan和Vese提出的M-S簡化模型，即C-V模型［8］。該隱性形狀先驗利用一個位置矢量和一個形狀矢量實現(xiàn)先驗形狀的變化。位置矢量用于空間變換，形狀矢量則是通過將符號距離函數(shù)隱性表示的形狀訓(xùn)練集進行主元分析（principal component analysis，PCA）得到且用于表示形狀差異的變化。2006年，Bresson等采用了和Tsai一樣的隱性先驗形狀表示方法［9］，不同之處在于，Bresson的模型在分割圖像的過程中，通過形狀、邊界和區(qū)域三個函數(shù)來調(diào)節(jié)圖像信息和形狀先驗知識在不同分割場合中的作用力度。

上述模型采用的先驗形狀模型或是固定的，或是可變的。不過即使是可變的，也只能在平均先驗形狀基礎(chǔ)上通過特征矢量的線性組合實現(xiàn)一些簡單的變化，變化范圍不大。因此，感興趣物體的形狀不能和先驗形狀相差太大。這是由于他們提取訓(xùn)練集特征的方法采用的是PCA。眾所周知，PCA是一種線性方法，只利用了先驗知識的二階統(tǒng)計信息。當(dāng)待分割物體的形狀過于復(fù)雜或者有非線性形變的時候，往往得不到滿意的分割結(jié)果。Cremers等2004年首先成功地將核方法通過Parzen估計應(yīng)用于幾何主動輪廓（geometric active contour，GAC）模型，從而構(gòu)建了形狀信息在核空間的分布模型［10］。Dambreville等2008年運用核主元分析（kernel PCA，KPCA）方法捕捉形狀信息，建立了形狀模型，并將其嵌入GAC分割模型中，成功地分割出正踢球的足球員和鯊魚等［11］。不同于PCA，KPCA能提取出數(shù)據(jù)的非線性特征，以及高階統(tǒng)計信息等更多的信息，應(yīng)用該模型能分割出形狀復(fù)雜的、或有非線性變形的目標(biāo)物體。

本研究提出了一種新的基于KPCA和形狀先驗知識的圖像分割模型。該模型中，分割曲線（含形狀訓(xùn)練集）由水平集隱式表示，利用KPCA提取訓(xùn)練集的非線性特征以擴大分割范圍；同時，將提取出的形狀先驗知識引入基于區(qū)域的可變形模型，運算效率高而且魯棒性強。這樣，能夠根據(jù)不同質(zhì)量的圖像，調(diào)節(jié)先驗知識和圖像信息兩者的作用力度，從而得到最佳的分割結(jié)果。

1 原理回顧

1.1 用水平集表示形狀

采用Leventon等2000年提出的隱式形狀表示法［7］，該方法是基于Osher和Sethian提出的水平集方法［2］。一個形狀/閉合曲線C通過一個更高維的符號距離函數(shù)φ隱式表示，φ是一種特殊的水平集函數(shù)，定義為

式中，d表示點（x，y）到曲線C的最短歐式距離，t是表示曲線進化過程中的時間步長。訓(xùn)練集中的形狀/曲線表示為高維函數(shù)φ的零水平截集，即C（t）=｛（x，y）|φ（x，y，t）=0｝。水平集表示法有很多優(yōu)點，如能精確地表示形狀，算法穩(wěn)定，易于實現(xiàn)曲線的分裂、合并、形成尖角等拓撲變化等。

1.2 KPCA

KPCA可以認為是線性PCA的推廣，是一種提取數(shù)據(jù)非線性特征的強有力法［12］。KPCA的基本思想很簡單，首先將輸入空間RN的數(shù)據(jù)X通過一個非線性變換φ映射到更高維的特征空間F，然后再在空間F中對這些投影數(shù)據(jù)進行主元分析（PCA）。一般來說，變換函數(shù)φ是非線性的，將本在輸入空間里不能線性分類的數(shù)據(jù)經(jīng)過映射變換后，在特征空間實現(xiàn)線性分類，如圖1所示。

圖1 線形空間RN到非線性空間F的映射［12］Fig.1 Projection from linearspaceRNto nonlinear oneF［12］

由于V屬于φ（X1），φ（X2），…，φ（Xl）的生成空間，因此有

并且存在參數(shù)αi，使得V可由φ（Xi）線性表出，即V=。定義一個n×n的矩陣K：Kij=k（Xi，Xj）=（φ（Xi）·φ（Xj）），將D、V的表達式代入式（2）中，于是可得式（1）的等價式為：nγα=Kα，α=（α1，α2，…，αn）T。因此求出K的特征值和特征向量，即求得式（1）的解。這樣就可提取出測試點X在映射空間的非線性主元，其第l主元為

2 基于KPCA的可變形模型

筆者所提出的是基于KPCA的可變形模型，該模型基于Chan和Zhu［13］提出的先驗形狀分割模型。設(shè)圖像矩陣I的領(lǐng)域范圍為Ω，距離函數(shù)φ表示進化曲線，ψ為訓(xùn)練集平均形狀的符號距離函數(shù)。源形狀φ和靶形狀ψ之間的關(guān)系用距離簡單地表示。在Chan-Vese（C-V）分割模型［8］和Paragios等［14］的形狀模型的基礎(chǔ)上，Chan-Zhu（C-Z）給出了其能量方程［13］

式中，H是Heviside函數(shù)，c1和c2分別為曲線C內(nèi)部和外部的圖像I的平均灰度值，λ≥0是圖像信息（第一項）和形狀先驗（第二項）的權(quán)重。形狀先驗項引導(dǎo)曲線φ逼近形狀先驗ψ，因此該模型能分割出部分遮擋的物體。

但是，先驗ψ是一個固定的形狀，為輸入空間的訓(xùn)練集的平均形狀，它僅僅吸引著曲線逼近平均形狀，這表明目標(biāo)物體的幾何形狀應(yīng)和先驗形狀的幾何形狀相接近。換言之，這是C-Z模型某種意義上的缺點，它沒有充分利用訓(xùn)練集里所包含的所有個體的差異先驗信息，而僅僅利用了他們的共有先驗信息——平均形狀，這極大地限制了該模型的分割范圍以及分割效果。

圖2 KPCA的逆映射Fig.2 The preimage problem of KPCA

采用的核函數(shù)為高斯核k（Xi，Xj）=exp（-‖Xi-Xj‖2/（2σ2）），那么圖像空間的距離函數(shù)=Eshape就可以用特征空間的距離函數(shù)來表示：

因此，基于KPCA的先驗形狀分割模型的能量方程為

然后運用歐拉方程和變分方法，對φ進行最小化E，可以得到曲線進化方程

參數(shù)λ和σ的選取關(guān)系到圖像分割質(zhì)量的高低。λ用于調(diào)節(jié)圖像信息和先驗信息在分割過程中分別所占的比重，σ的大小影響著所獲得的曲線形狀和先驗形狀的差異程度。因此選取時，要綜合考慮圖像的質(zhì)量和物體的幾何形狀。當(dāng)圖像質(zhì)量不好時（比如物體部分被遮擋或者有殘缺，或者噪聲太大，或者有其他灰度值接近的物體和感興趣物體相連），應(yīng)該由先驗信息來正確引導(dǎo)曲線的進化。因此應(yīng)增大λ的值以加強先驗項的比重。同時應(yīng)該減小σ，將曲線的形狀嚴格的限制在一定的形變范圍內(nèi)，減少其他因素的影響。反過來，當(dāng)感興趣目標(biāo)的幾何形狀和先驗形狀相差較大時，要得到合適的有意義的進化曲線，應(yīng)更多的考慮給定的圖像信息，而不是先驗信息。因此應(yīng)減小λ增強圖像項的比重，同時增大σ擴大曲線的形變范圍。此時曲線的進化結(jié)果更多的依賴于給定的圖像信息。

3 算法驗證

利用合成圖像和真實圖像，將筆者的模型和基于PCA的分割模型分割結(jié)果進行比較，對所提出的算法進行驗證。

3.1 采用合成圖像

采用文獻［6］中的飛機訓(xùn)練集。該訓(xùn)練集包含12架飛機（如圖3（a）所示），每幅圖像大小為231像素×231像素。由于飛機的大小位置相互間存在差異，在提取形狀先驗前需進行配準(zhǔn)，從而消除大小位置帶來的誤差。配準(zhǔn)方法也采用文獻［6］中的方法，第一架飛機固定不變，作為參考，將其余飛機進行剛性變換，與其對齊，得到圖3（b）所示的配準(zhǔn)結(jié)果。

另外，選1幅不包含于訓(xùn)練集的飛機圖像用于測試實驗，并在這幅測試圖的基礎(chǔ)上，再添加了其他的干擾因素，則生成3幅待分割圖像。第1幅中的飛機被大量的高斯噪聲所污染，噪聲信噪比為10dB；第2幅為有殘缺的飛機，其50～60行被置為背景色，并添加了10dB的高斯噪聲；第3幅采用的干擾因素為3×5柵欄和10dB高斯噪聲。每根條柵的寬度為7個像素，灰度值與飛機一致。豎條柵的起始橫坐標(biāo)依次為35、75、115、150和190；橫條柵的起始縱坐標(biāo)依次為40、105和180。在柵欄的干擾下，飛機被嚴重污染，且部分邊界被模糊。

考慮到基于PCA的分割模型的初始曲線不能和分割結(jié)果相差太大，一般選用平均形狀作為初始曲線。因此在下面的分割實驗中所選用的初始曲線均為平均形狀。

3.2 真實圖像

選取36張大小為128像素×128像素的含肝CT圖像，并由專家手工分割出肝，將此作為提取肝形狀特征的訓(xùn)練集。

同樣地，用于分割測試的肝圖像不包含于訓(xùn)練集，初始曲線為平均形狀。

4 結(jié)果

4.1 合成圖像

圖3中（a）和（b）分別為采用文獻［6］中的飛機訓(xùn)練集和相應(yīng)的配準(zhǔn)結(jié)果?？梢钥闯?，經(jīng)過配準(zhǔn)后，圖（b）中的飛機的大小和位置基本一致，機頭均朝正上方。

圖3 含有12架飛機的訓(xùn)練集。（a）配準(zhǔn)前；（b）配準(zhǔn)后Fig.3 Training data of 12 binary airplanes.（a）before alignment；（b）after alignment

圖4～圖6分別為人工合成噪聲污染圖像、殘缺圖像和柵欄污染圖像及其分割結(jié)果。各圖中，（a）為包含初始曲線的測試圖，（b）和（c）為初始曲線分別在PCA分割模型和本模型驅(qū)使下演變的最后結(jié)果?？傮w而言，本模型在3幅圖像上得到的結(jié)果均好于基于PCA分割模型的結(jié)果。盡管基于PCA分割模型幾乎檢測出了整個飛機，但是仍有部分曲線和飛機的實際邊界有一定的距離，如飛機頭，兩側(cè)的翅膀和機尾。特別是圖6中分割出的翅尖不但包含真實的邊界，還包含與其相連的背景柵欄。而本模型則給出了令人滿意的結(jié)果，特別是圖4，通過調(diào)節(jié)圖像信息項和先驗信息項在分割中的比重，使得曲線最終吸引到正確的物體邊界。即使是存在強干擾的圖5和圖6，本模型也給出了不錯的結(jié)果，檢測出了想要的飛機形狀。這得歸功于KPCA，一方面KPCA將分割曲線限定在先驗形狀周圍，排除外界的干擾，另一方面又允許分割曲線發(fā)生較大的或非線性的形變，更貼近飛機真實的邊界。

圖4 噪聲污染圖像的分割結(jié)果。（a）初始曲線；（b）PCA分割模型；（c）本模型（λ=0.9，σ=1）Fig.4 Segmentation of the noised airplane.（a）initialization；（b）result using PCA；（c）result using our model

圖5 殘缺圖像的分割結(jié)果。（a）初始曲線；（b）PCA分割模型；（c）本模型（λ=0.8，σ=3）Fig.5 Segmentation of the incomplete airplane.（a）initialization；（b）result using PCA；（c）result using our model

4.2 真實圖像

與合成圖像同理，需在提取形狀先驗前將36幅肝CT圖像訓(xùn)練集進行配準(zhǔn)。因篇幅原因，圖7僅列出了其中的15幅，如圖7所示。其中，（a）為配準(zhǔn)前的圖像集，（b）為配準(zhǔn)結(jié)果，圖中肝的大小、位置和朝向幾乎一致。

圖6 柵欄污染圖像的分割結(jié)果。（a）初始曲線；（b）PCA分割模型；（c）本模型（λ=0.8，σ=3）Fig.6 Segmentation of the grilled airplane.（a）initialization；（b）result using PCA；（c）result using our model

圖8所示為CT肝圖像的分割結(jié)果?？梢钥闯?，PCA分割模型分割出了肝的整體形狀，不過在某些肝的細節(jié)上卻差強人意，如右上角的尖角和中下方的凹處都沒有得到很好的分割。相比較而言，本模型在分割肝的過程中，一方面KPCA允許進化曲線的形變范圍擴大，另一方面肝的灰度信息和邊界信息正確地引導(dǎo)進化曲線逼近正確的肝的邊界，因此最終在這兩處都得到了準(zhǔn)確的分割結(jié)果。

圖7 含有15個肝的訓(xùn)練集。（a）配準(zhǔn)前；（b）配準(zhǔn)后Fig.7 Training data of 15 binary livers.（a）before alignment；（b）after alignment

圖8 CT肝圖像的分割結(jié)果。（a）初始曲線；（b）PCA分割模型；（c）本模型（λ=0.8，σ=1）Fig.8 Segmentation of the liver.（a）initialization；（b）result using PCA；（c）result using our model

5 討論

由于成像過程中存在各種噪聲污染、電磁等影響，醫(yī)學(xué)圖像質(zhì)量不高，普遍存在邊緣模糊、殘缺、不同組織之間的界線不清晰等問題。而且呼吸運動和個體差異也導(dǎo)致欲分割的器官存在一定的形變和差異。此時僅利用圖像信息的可形變模型很難得到理想的分割結(jié)果。如圖9（c）所示，僅利用圖像信息的C-V模型盡管分割出了肝，但也分割出了和肝灰度值近似的其他器官如胃、脾等。

對此，筆者所提出的模型融入了用KPCA分析的形狀先驗知識，利用KPCA代替PCA能提出更有價值的先驗形狀的特征向量，擴大了分割范圍，使得分割模型不但可以正確處理一般的污染嚴重、有遮擋和信息缺失的情況；而且，即使待分割形狀與先驗形狀有著非線性的形變，也可以得到正確分割。在KPCA分析得到的形狀模型引導(dǎo)下，分割質(zhì)量不好、存在個體差異、卻有著豐富訓(xùn)練集的醫(yī)學(xué)圖像，也能得到滿意的分割結(jié)果。在圖8的分割結(jié)果中，盡管平均形狀的右下方?jīng)]有凹處，但KPCA允許進化曲線發(fā)生較大的形變，因此最終的結(jié)果出現(xiàn)了本研究想要的凹處，而PCA則沒有出現(xiàn)。不過KPCA的該性能也使得右上角分割出了和肝緊連的部分胃壁，這和本研究采用的C-V模型有關(guān)，該模型所采用的分割特征為灰度值，因此灰度值接近的胃壁也被分割出來了。在分割過程中，雖然可以調(diào)節(jié)σ和λ來得到不同的分割結(jié)果，但想得到很滿意的分割結(jié)果，參數(shù)的調(diào)節(jié)并不容易。

圖9 CT肝圖像的分割結(jié)果。（a）原始圖像；（b）專家分割結(jié)果；（c）C-V模型Fig.9 Segmentation of the liver.（a）original image；（b）result by expert；（c）result using C-V model

再則，該模型選取了任意的一個逆映射函數(shù)，而不是固定的平均形狀，作為參考形狀。該逆映射函數(shù)是分割進化曲線在KPCA特征子空間的映射，它隨著曲線的進化而發(fā)生著變化，最終的進化曲線可以和訓(xùn)練集中任意一個個體近似，充分利用了訓(xùn)練集里所包含的所有個體的差異先驗信息。這給初始曲線的選擇預(yù)留了足夠的空間，不僅僅局限于平均形狀。圖10中，利用手工繪制的多邊形作為初始曲線，盡管該曲線和最終曲線的形狀相差較大，但由于進化過程中該初始曲線不斷地逼近訓(xùn)練集中的個體，因此得到的最終結(jié)果還是令人滿意的。

圖10 選用其他初始曲線的分割結(jié)果。（a）初始曲線；（b）分割結(jié)果（λ=0.7，σ=2）Fig.10 Result using another initialization.（a）initialization；（b）result

6 結(jié)論

僅利用圖像信息的可形變模型對于質(zhì)量不高的醫(yī)學(xué)圖像來說很難得到理想的分割結(jié)果。融合先驗知識和圖像信息是提高分割質(zhì)量的一條有效途徑。本研究提出了一種新的基于形狀先驗知識和KPCA的圖像分割模型，并將該模型應(yīng)用于醫(yī)學(xué)CT圖像。該模型利用KPCA提取先驗知識，構(gòu)建形狀模型，并將其融入基于區(qū)域的C-V模型。在比較實驗中，所提出的分割模型應(yīng)用于合成圖像和醫(yī)學(xué)CT圖像，取得了滿意的結(jié)果。與基于PCA的模型所得的分割結(jié)果相比較，該模型更能準(zhǔn)確地識別待分割物體，尤其當(dāng)圖像質(zhì)量較差（背景污染嚴重、缺損、含錯誤信息等）和待分割物體的形狀與先驗形狀差別較大（非線性形變等）時具有顯著優(yōu)勢。

本研究所提出的分割模型在進行KPCA分析時，采用了經(jīng)典的高斯核函數(shù)，其他核函數(shù)（如多項式核函數(shù)）在應(yīng)用于圖像分割中的性能表現(xiàn)還值得進一步研究。此外，本研究的目標(biāo)是將此模型進一步擴展應(yīng)用到3D醫(yī)學(xué)圖像，使其自動分割更具準(zhǔn)確性和魯棒性。

［1］Kass M，Witkin A，Terzopoulos D.Snakes：Active contour models［J］.International Journal of Computer Vision，1988，1（4）：321-331.

［2］Osher S，Sethian J.Fronts propagating with curvature dependent speed：algorithms based on the hamilton-jacobi［J］.Journal of Computational Physics，1988，79（1）：12-49.

［3］Mumford D，Shah J.Optimal approximations by piecewise smooth functions and associated variational problems［J］.Communications on Pure and Applied Mathematics，1989，42（5）：577-685.

［4］Chen Y，Thiruvenkadam S，Huang F，etal.On the incorporation of shape priors into geometric active contours［C］//Proceedings of the IEEE Workshop on Variational and Level Set Methods.Washington：IEEE，2001：145-152.

［5］Cremers D，Tischhauser F，Weickert J，et al.Diffusion snakes：Introducing statistical shape knowledge into the mumford-shah functional［J］.International Journal of Computer Vision，2002，50（3）：295-313.

［6］Tsai A，Yezzi A，Wells，W，et al.A shape-based approach to the segmentation of medical imagery using level sets［J］.IEEE transactions on medical imaging，2003，22（2）：137-154.

［7］Leventon M，Grimson W，F(xiàn)augeras O.Statistical shape influence in geodesic active contours［C］//Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition.South Carolina：IEEE，2000：316-323.

［8］Chan T，Vese L.Active contours without edges［J］.IEEE Transactions on Image Processing.2001，10（2）：266-277.

［9］Bresson X，Vandergheynst P，Thiran J.A variational model for object segmentation using boundary information and shape prior driven by the mumford-shah functional［J］.International Journal of Computer Vision，2006，68（2）：145-162.

［10］Cremers D，Osher S，Soatto S.Kernel density estimation and intrinsic alignment for knowledge-driven segmentation：teaching level sets to walk［C］//Proceedings of Springer Pattern Recognition.Rasmussen：Springer，2004：36-44.

［11］Dambreville S，Rathi Y，Tannenbaum A.A framework for image segmentation using shape models and Kernel space shape priors［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2008，30（8）：1385-99.

［12］Sch?lkopf B，Smola A，Müller K.Nonlinear component analysis as a kernel eigenvalue problem［J］.Neural Computation，1998，10（5）：1299-1319.

［13］Chan T，Zhu W.Level set based shape prior segmentation［C］//Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition.Washington，DC：IEEE，2005：1164-1170.

［14］ParagiosN，Rousson M，Ramesh V.Matching distance functions：A shape-to-area variational approach for global-tolocal registration［C］//Proceedings of Computer Vision.Copenhagen：Springer，2002：813-815.

［15］Kwok J，Tsang I.The pre-image problem in kernel methods［J］.IEEE Trans Neural Netw，2004，15（6）：1517-1525.