基于瓶頸檢測(cè)和分水嶺算法的重疊宮頸細(xì)胞圖像分割方法

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

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

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

0 引言

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

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

1 方法

1.1 預(yù)處理

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

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

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

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

1.2 確定分裂點(diǎn)

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

1.2.1 多邊形近似

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

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

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

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

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

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

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

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

圖2 檢測(cè)特征點(diǎn)Fig.2 Feature detection

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

1.2.2 瓶頸檢測(cè)

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

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

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

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

圖3 瓶頸檢測(cè)的例子Fig.3 Examples of bottleneck detection

1.2.3 橢圓擬合

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

首先，把所有候選分裂點(diǎn)對(duì)根據(jù)瓶頸率從低到高排序。從第一對(duì)候選分裂點(diǎn)開(kāi)始，將重疊細(xì)胞區(qū)域的輪廓?jiǎng)澐譃長(zhǎng)1和L2兩段。然后，利用多邊形近似得到的每個(gè)線段上的特征點(diǎn)擬合橢圓，表示為：

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

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

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

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

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

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

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

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

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

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

圖6 分割過(guò)程Fig.6 The segmentation process

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

2 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

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

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

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

3 討論和結(jié)論

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