段鵬，程文播，錢慶，章強(qiáng)，楊任兵，潘宇駿

1 中國科學(xué)技術(shù)大學(xué)，合肥市，230026

2 中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所，蘇州市，215163

0 引言

國家癌癥中心2018年統(tǒng)計的數(shù)據(jù)表明，在中國，每年被檢測出的宮頸癌患者達(dá)到13萬，因?qū)m頸癌去世的女性人數(shù)約5.2萬人。宮頸癌是危害女性健康和生命的惡性腫瘤之一。宮頸癌的篩查工作對廣大婦女的生命健康有著非凡的作用[1]。目前經(jīng)常被使用的宮頸癌的篩查方法是宮頸細(xì)胞學(xué)檢測[2]?？紤]到宮頸癌普查的工作量極大，通過計算機(jī)圖像處理技術(shù)輔助宮頸癌的普查具有重要的現(xiàn)實意義[3]。

宮頸細(xì)胞的形態(tài)、大小等特征是用來判斷宮頸疾病的指標(biāo)之一[4]。宮頸細(xì)胞圖像中，細(xì)胞質(zhì)之間往往會發(fā)生大面積重疊，從重疊的細(xì)胞區(qū)域分割出單獨的宮頸細(xì)胞區(qū)域是一項極具挑戰(zhàn)性的任務(wù)[5]，染色過程中引起的光照差異和染色濃度等問題都會導(dǎo)致細(xì)胞不同成分之間的顏色對比度降低。其它一些例如空氣干燥、過多的血液、細(xì)菌、炎癥的問題也會使得重疊細(xì)胞區(qū)域的識別更具挑戰(zhàn)性[6]。此外，宮頸細(xì)胞還存在本身的形狀相對比較復(fù)雜，重疊程度較高，細(xì)胞質(zhì)的邊界信息較弱等問題。因此如何在宮頸細(xì)胞的篩查過程中，提高宮頸細(xì)胞圖像的分割效果，將為宮頸癌疾病的治療提供了最有效的幫助。因此本次實驗的目的在于討論重疊宮頸細(xì)胞圖像的分割算法。

1 方法

1.1 預(yù)處理

本研究使用最大類間方差算法對重疊宮頸細(xì)胞圖像進(jìn)行預(yù)處理，獲取宮頸細(xì)胞的輪廓信息。最大類間方差算法簡稱Otsu算法，是一種自動適應(yīng)確定閾值的方法。算法實現(xiàn)假設(shè)圖像的所有像素一定可以根據(jù)閾值被準(zhǔn)確分成細(xì)胞區(qū)域和背景區(qū)域兩部分。在此基礎(chǔ)上，若需要得到更準(zhǔn)確的分割結(jié)果，那么就要盡可能使得細(xì)胞區(qū)域與背景區(qū)域這兩部分的灰度值區(qū)分度越大，而方差正好是用來描述灰度值分布是否均勻的方法，因此通過計算某個閾值處，兩個區(qū)域的類間方差，當(dāng)類間方差值達(dá)到最大時，此時的閾值就是能夠使圖像最準(zhǔn)確分割的閾值點。

假設(shè)一個圖像I(x，y)，它的像素點區(qū)域為M×N，將閾值記為T，假定此時細(xì)胞區(qū)域包含的所有像素的總和占所有像素的比重是w0，像素的平均值是μ0，同理，假定背景區(qū)域包含的所有像素的總和占所有像素的比重為w1，像素的平均值是μ1。整個圖像灰度值的平均值設(shè)為μ，方差設(shè)為g。此外，我們將灰度值低于閾值T的像素總和設(shè)為N0，灰度值高于閾值T的像素總和設(shè)為N1，有如下的關(guān)系：

將式(5)代入式(6)，得到：

將閾值從最小依次調(diào)整到最大，每次設(shè)定閾值后，求出類間方差g，當(dāng)g的值最大時，此時預(yù)先設(shè)定閾值T就是使圖像完美分割的最佳預(yù)置點，使用該閾值可以將細(xì)胞區(qū)域和背景區(qū)域分割出來。

1.2 確定分裂點

重疊細(xì)胞在重合的地方往往都會出現(xiàn)凹陷的現(xiàn)象，兩個細(xì)胞之間會出現(xiàn)凹點而且是成對出現(xiàn)，因此找到成對的凹點是進(jìn)行重疊細(xì)胞正確分割的關(guān)鍵點[7]。我們將這個凹點對稱之為使重疊細(xì)胞分割的分裂點。大多數(shù)傳統(tǒng)的凹點檢測方法是基于凸圖形的定義提出的一種簡便的凹點檢測方法，原理簡單，但需要大量冗余的計算，耗時較長?？紤]到想讓算法盡可能快速高效地處理圖像，本研究結(jié)合多邊形近似、瓶頸檢測和橢圓擬合進(jìn)行分裂點檢測。

1.2.1 多邊形近似

多邊形近似[8]的目的是在某一曲折連續(xù)的線段上進(jìn)行判斷，根據(jù)多邊形的特點檢測出特征點作為該線段上近似于多邊形的頂點，從而找出可能存在重疊細(xì)胞的特征點區(qū)域。一個理想的多邊形近似的算法應(yīng)當(dāng)以盡可能少的的計算復(fù)雜度，卻盡可能地標(biāo)記處出線段的所有形狀特點。本研究采用的多邊形近似的方法如下：

假設(shè)C={p1,p2,...,pi,...,pn}代表曲線的原始輪廓，式中，pi=(xi,yi)代表曲線上的每一個點，n代表原始輪廓上的點數(shù)。

第一步：隨機(jī)選擇兩個輪廓點pi(i=1)和pj(j=i+2)，如圖1所示：

圖1 多邊形近似的例證Fig.1 Example of polygon approximation

第二步：計算pi和pj之間每一個輪廓點到lpi pj的距離，并和預(yù)設(shè)閾值Td進(jìn)行比較，如圖1所示。如果存在點pt，其對應(yīng)的距離dt滿足dt=Td，那么pt被認(rèn)為是輪廓的分裂點，并且pi和pj從pt和pt+2開始繼續(xù)，否則，pi保持不變，而pj變成pj+1繼續(xù)。

第三步：重復(fù)第二步，直到遍歷了整個輪廓的所有點。

在多邊形近似中，閾值Td用于控制精度，閾值設(shè)置越小越好。但是相應(yīng)帶來的缺點就是：當(dāng)閾值設(shè)置過小，那么曲線上得到的特征點就過多。因此，在后續(xù)步驟中的計算時間也會變的越長。所以，閾值Td出于效率和精度兩方面的考慮。

該多邊形近似的效果如圖2所示，圖左為某一細(xì)胞的輪廓信息，圖右是使用上述步驟得到的所有特征點標(biāo)注。

圖2 檢測特征點Fig.2 Feature detection

多邊形近似處理后得到的所有特征點，都有可能是兩個細(xì)胞重合部分的分裂點。因此，需要從所有特征點中找到正確的分裂點，才能正確地分割重疊細(xì)胞。

1.2.2 瓶頸檢測

瓶頸檢測[9]的主要目的是從上述多邊形近似找到的所有特征點中，找到兩個重疊細(xì)胞重合的凹點，即能夠使重疊細(xì)胞被正確分割的分裂點，正確識別重疊細(xì)胞的分裂點是實現(xiàn)重疊細(xì)胞正確分割的重要步驟。

瓶頸檢測的方法大多數(shù)都是基于凹面和形狀的算法，本文瓶頸檢測算法如下所示：對于通過多邊形近似得到的所有特征點，通過計算每對特征點之間的瓶頸率來判斷該特征點對是否為正確的分裂點對。假設(shè)兩個特征點pa和pb，瓶頸率e的計算方法是：

式中，dist(pa,pb)為pa和pb兩點間的歐幾里得距離，length(pa,pb)為點pa到pb順時針方向的邊界長度。length(pa,pb)和length(pb,pa)分別表示pa和pb兩點之間的上邊緣和下邊緣的長度。若pa和pb兩個特征點計算出的瓶頸率小于閾值Te，那么這兩個特征點就被認(rèn)定為是一組候選分裂點對。對于標(biāo)準(zhǔn)圓形單元，瓶頸率的最小值為2/π。但是，考慮到宮頸細(xì)胞圖像可能存在不規(guī)則性或橢圓形，閾值Te應(yīng)該適當(dāng)小于2/π，因此本文的閾值設(shè)為1/2。如果沒有檢測到候選分裂點對，那表明該區(qū)域只包含一個細(xì)胞單元。

該算法可以找到一些候選分裂點對，但當(dāng)單細(xì)胞圖像本身呈現(xiàn)外凸時，也有可能檢測到分裂點對，而此時的分裂點對是錯誤的分裂點對。圖3就是一個未正確檢測到分裂點對的例子，右圖的細(xì)胞輪廓為一個單細(xì)胞，此時點與點之間的瓶頸率為0.30，滿足e(pa,pb)＜Te，算法判斷它是存在分裂點對，出現(xiàn)了明顯的誤判斷。這就說明利用瓶頸檢測檢測出來的分裂點對不一定是正確的分裂點對。

圖3 瓶頸檢測的例子Fig.3 Examples of bottleneck detection

1.2.3 橢圓擬合

橢圓擬合[10]主要有聚類和最小二乘擬合兩大類方法。相比聚類方法，最小二乘法的橢圓擬合技術(shù)擁有耗時更少、魯棒性更好、精度更高等優(yōu)勢。鑒于此，本研究采用一種被廣泛應(yīng)用的最小二乘擬合方法直接作為橢圓擬合的工具。

首先，把所有候選分裂點對根據(jù)瓶頸率從低到高排序。從第一對候選分裂點開始，將重疊細(xì)胞區(qū)域的輪廓劃分為L1和L2兩段。然后，利用多邊形近似得到的每個線段上的特征點擬合橢圓，表示為：

通常，宮頸細(xì)胞在圖像中的區(qū)域范圍差異應(yīng)該不大。因此，如果兩個擬合橢圓的面積差較大時，那么該分裂點對可能是不正確的。用S1和S2分別表示根據(jù)L1和L2上的特征點擬合的兩個橢圓的面積，r是由r=max(S1/S2,S2/S1)計算出的面積比的最大值，如果滿足r＞Tr，則認(rèn)為此時對應(yīng)的候選分裂點不正確。其中Tr是一個預(yù)設(shè)閾值，另外，還利用平均代數(shù)距離（MAD）作為另一種度量，進(jìn)一步確定候選分裂點對是否正確。令xi代表線段L上的特征點(xi,yi),i=1,2,...,s，在線段L上所有特征點與擬合橢圓之間的MAD計算公式[11]如下：

式中，s代表線段L上所有特征點的數(shù)目，xi=[xi2,代表點(xi yi)到橢圓ax2+bxy+cy2+dx+ey+f=0的距離，根據(jù)+f計算得出。

分別計算線段L1和L2的MAD1和MAD2，并將它們的最小值MADmin與預(yù)先設(shè)定的閾值Tmad進(jìn)行比較，如果MADmin＜Tmad，表示至少有一個分段L1或者L2的擬合誤差要小得多，此時可以認(rèn)為對應(yīng)的候選分裂點對是正確的。對于偏心率較大的單細(xì)胞，瓶頸檢測可能會誤檢出一些分裂點對。此時在這種情況下，這兩個條件的擬合誤差都比較小，因此不能用MADmin＜Tmad消除不正確的分裂點對。為了解決這一問題，增加了另一項限制為：MADc＜Tmad，其中，MADc表示從整個輪廓計算出的平均代數(shù)距離。從單個細(xì)胞中計算出的MADc值通常會比較小，因為它們近似于橢圓形狀。單個細(xì)胞可以通過閾值Tmad來比對它們的擬合誤差MADc，從而去除不正確的分裂點對。只有當(dāng)MADmin和MADc分別滿足上述的限制條件時，才能將對應(yīng)區(qū)域視為重疊細(xì)胞區(qū)域，并判斷此時的候選分裂點對是正確的分裂點對。

圖4中顯示了錯誤分割的例子和使用橢圓擬合校正分裂點對，其中實點代表細(xì)胞圖像輪廓上的特征點，虛線代表使用特征點擬合成的橢圓形狀，其中特征點pa和pb表示候選分裂點對。在圖4的左半部分，擬合橢圓對特征點的擬合不正確，擬合橢圓計算得到的MADmin值為0.031，即該細(xì)胞輪廓對應(yīng)的區(qū)域不是一個單細(xì)胞，認(rèn)為此時的候選分裂點對不正確，應(yīng)予以丟棄。作為對比，給出了一個用正確的分裂點對擬合橢圓的例子，其中MADc=0.171，MADmin=0.009?？梢钥闯觯卣鼽c與擬合橢圓的分布非常吻合，擬合橢圓可以幫助找到正確的分裂點對，從而能夠正確分割相應(yīng)的區(qū)域。如果當(dāng)前的候選分裂點對判斷為不正確，則將丟棄該候選分裂點對，并選擇下一對進(jìn)行判斷，直到找到一個滿足條件的分裂點對為止。如果所有的分裂點對都不滿足條件，則表示該重疊區(qū)域其實是單個細(xì)胞區(qū)域。

圖4 橢圓擬合校正Fig.4 Correction by ellipse fitting

1.3 重疊細(xì)胞重合區(qū)域的分割

根據(jù)1.2節(jié)的算法，可以得到重疊細(xì)胞的正確分裂點對，由于宮頸細(xì)胞大多數(shù)呈現(xiàn)圓形或者似圓形。因此，我們可以斷定重疊細(xì)胞的重合區(qū)域應(yīng)當(dāng)包含于以兩個分裂點連接的線段為直徑的圓內(nèi)，如圖5所示。

因為在梯度圖像中，細(xì)胞邊緣處的灰度值通常較高，并且能夠很好地保留下細(xì)胞圖像。1.2節(jié)中得到了重疊細(xì)胞的正確分裂點對，根據(jù)分裂點對又可以縮小重疊細(xì)胞重合的區(qū)域，優(yōu)化算法難度。本研究先求出兩個重疊細(xì)胞的重合部分的梯度圖像，再采用分水嶺算法對該重合區(qū)域的梯度圖像進(jìn)行分割，這樣可以大大縮小圖像處理的范圍，從而可以提高分割的性能。但是，直接用分水嶺算法處理梯度圖像往往會出現(xiàn)過分割的情況，因此，本研究預(yù)先使用閾值算法對梯度圖像進(jìn)行處理，去除掉微弱（梯度較?。┑倪吔缧畔ⅲD像的梯度圖像歸一化到區(qū)域[0，1.0]之間，再對梯度圖像進(jìn)行分水嶺算法，從而可以求出重疊細(xì)胞重合區(qū)域的分割邊界，再與預(yù)處理獲得的細(xì)胞前景區(qū)域輪廓進(jìn)行疊加，從而最后得到完整的重疊細(xì)胞的分割效果。

圖5 重合區(qū)域的確定Fig.5 Determination of coincidence region

由圖6可以看出，該算法明顯地去除了重疊細(xì)胞重合區(qū)域的微弱邊界區(qū)域，同時為了解決過分割的問題，將重疊細(xì)胞的重合區(qū)域單獨拿出來進(jìn)行分割，得到重合區(qū)域的邊界后，再疊加用細(xì)胞前景分割算法提取出的外細(xì)胞輪廓，構(gòu)成了完整的重疊細(xì)胞的分割結(jié)果。

圖6 分割過程Fig.6 The segmentation process

綜上所述，重疊細(xì)胞分割算法步驟如下：

2 實驗結(jié)果與分析

本實驗使用華碩臺式電腦（Intel i5-4590 3.3 GHz CPU處理器，8 GB內(nèi)存，GTX1050Ti顯卡），使用Matlab R2017a軟件作為開發(fā)工具進(jìn)行研究開發(fā)。

實驗部分將50幅重疊的宮頸細(xì)胞顯微圖像進(jìn)行分割，圖像尺寸為250×250，為了說明算法的分割結(jié)果，將分割結(jié)果與醫(yī)生手工分割結(jié)果進(jìn)行比較，并與文獻(xiàn)[9]的分割結(jié)果進(jìn)行對比。部分分割結(jié)果如圖7所示，其中第一行為原始輸入圖像，第二行為本研究算法的分割結(jié)果，第三行為醫(yī)生標(biāo)注的重疊細(xì)胞的分割結(jié)果。

圖7 重疊細(xì)胞分割結(jié)果Fig.7 Overlapping cell segmentation results

從圖7的分割結(jié)果可以得出，本研究算法能夠?qū)⒅丿B的宮頸細(xì)胞圖像有效分離開來，實現(xiàn)了高重疊細(xì)胞的顯微圖像分割，且分割效果接近于醫(yī)生的手工分割結(jié)果。針對重疊區(qū)域輪廓較為明顯的重疊細(xì)胞，本研究的分割算法較好，前三列的分割結(jié)果與醫(yī)生手工分割結(jié)果接近，誤差較小。最后一列的重疊區(qū)域顏色較深，輪廓較為模糊，難以分割準(zhǔn)確。但相對文獻(xiàn)[9]已經(jīng)有了一定的提高。本研究采用分割區(qū)域面積的誤差比w來定量評價分割結(jié)果，算式為：

式中，areas為本研究分割出的區(qū)域面積，arear為醫(yī)生手工分割的區(qū)域面積。

本研究選取目前國內(nèi)外主流的宮頸細(xì)胞圖像分割算法進(jìn)行實驗，并將他們的結(jié)果進(jìn)行對比分析，如表1所示。

表1 本研究算法與其它算法的比較Tab.1 Comparison of proposed algorithm and other algorithms

由表1可知，雖然本研究算法運(yùn)行時間較長，犧牲了一定的運(yùn)行效率。但本研究算法的分割結(jié)果的面積誤差為0.08，是四種算法中效果最好的，說明與醫(yī)生的手工分割結(jié)果較為接近。

3 討論和結(jié)論

針對重疊的宮頸細(xì)胞圖像重疊情況復(fù)雜、難以分割的問題。本研究提出了一種基于瓶頸檢測和分水嶺算法的分割算法。首先，通過細(xì)胞前景區(qū)域的分割算法獲得整個重疊細(xì)胞的輪廓。其次，通過瓶頸檢測和橢圓擬合來檢測重疊細(xì)胞輪廓的分裂點，從分裂點中判斷出正確分割圖像的分裂點對。再次，根據(jù)分裂點進(jìn)一步縮小重疊細(xì)胞重合的區(qū)域范圍，降低處理難度。最終通過分水嶺算法獲得重疊區(qū)域的邊界信息。將重疊區(qū)域的邊界與外輪廓疊加起來，獲得重疊細(xì)胞的分割結(jié)果。將該算法應(yīng)用于重疊的宮頸細(xì)胞圖像分割，能夠分割出準(zhǔn)確完整的單個宮頸細(xì)胞圖像，且與醫(yī)生手動標(biāo)注的分割結(jié)果接近，與最新的其他基于瓶頸檢測的分割算法相比，誤差更小，更為準(zhǔn)確。但是犧牲了一定的運(yùn)行效率，且當(dāng)重疊區(qū)域非常復(fù)雜時，仍然存在較大的分割誤差，如何解決這兩個問題有待進(jìn)一步的研究改進(jìn)。