陳書文,曹 愚,趙小燕,王茄吉

(1.江蘇第二師范學(xué)院 數(shù)學(xué)與信息技術(shù)學(xué)院，江蘇 南京 210013；2.南京工程學(xué)院 信息與通信工程學(xué)院，江蘇 南京 211167)

細(xì)胞數(shù)字圖像是臨床醫(yī)學(xué)、細(xì)胞學(xué)、病理學(xué)的重要研究手段，在疾病診斷、癌細(xì)胞篩查中發(fā)揮著重要作用。憑借計(jì)算機(jī)圖像分析技術(shù)，醫(yī)生或科研人員能夠?qū)δ繕?biāo)細(xì)胞進(jìn)行定性或定量分析。這不僅能減少視覺工作量，還極大地提高了工作效率。細(xì)胞圖像分割是細(xì)胞識別與統(tǒng)計(jì)的核心，當(dāng)前分割方法可大致歸納為四類，第一類為閾值法；第二類為基于區(qū)域的方法；第三類為基于先驗(yàn)?zāi)Ｐ偷姆椒?；第四類是基于神?jīng)網(wǎng)絡(luò)的方法。

閾值法[1]把每個(gè)像素的灰度值作為特征，通過比較灰度值對像素進(jìn)行分類，一般分為單閾值法和多閾值法[2-3]，具有簡單易行、性能穩(wěn)定等優(yōu)點(diǎn)[4]。目前流行的閾值法有Otsu法[5]、最大熵法[6]和聚類法[7]等。此類方法對于簡單或無噪聲的細(xì)胞圖像表現(xiàn)較好，若遇光照不均勻、染色不均勻、或污染干擾的情況，細(xì)胞分割的準(zhǔn)確性將受到影響。基于區(qū)域的方法不僅要考慮圖像的灰度值，還要考慮融合環(huán)境。此類方法有區(qū)域生長法[8-9]、分水嶺法[10-11]等。區(qū)域生長法要先確定生長準(zhǔn)則，而分水嶺法的缺點(diǎn)是容易過分割，這都不能獨(dú)立地用于分析形狀小、密度高的細(xì)胞圖像。基于先驗(yàn)?zāi)Ｐ偷姆椒ㄊ峭ㄟ^預(yù)定的形狀來提取特定細(xì)胞，再進(jìn)行分類。但實(shí)際問題是細(xì)胞形態(tài)各異，預(yù)先確定細(xì)胞模板很困難。而基于神經(jīng)網(wǎng)絡(luò)的方法是把圖像的分割作為函數(shù)的最小化問題來處理，主要思想是用已知的結(jié)果作為樣本對神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練，此方法過程也較復(fù)雜。為尋找快速有效、且能應(yīng)對各種染色污染問題的細(xì)胞圖像分割法，本文提出基于K-means聚類和Canny算子相結(jié)合的新思路，并把新算法的效果與單獨(dú)使用以上兩種方法進(jìn)行對比，旨在證明新方法的優(yōu)越性。

1 算法原理

圖1 算法流程圖Fig.1 Flow diagram of algorithm

給定的細(xì)胞圖像經(jīng)過灰度化和中值濾波后，采用K-means聚類作為主分割法，同時(shí)用Canny算子的分割結(jié)果作為必要補(bǔ)充，把因染色問題導(dǎo)致未被正確檢測的細(xì)胞統(tǒng)計(jì)進(jìn)來，發(fā)揮兩種算法的優(yōu)點(diǎn)，提高細(xì)胞識別的準(zhǔn)確率(算法流程如圖1所示)。要注意的是，若原圖像為RGB彩色圖像，需先灰度化。預(yù)處理時(shí)提取圖像的RGB三個(gè)通道值，選擇直方圖均衡最好的通道建立灰度圖像。一般情況，圖像中的細(xì)胞核、細(xì)胞質(zhì)、環(huán)境背景是3種不同的灰度值，所以首次分割采用的K-means聚類，聚類中心數(shù)設(shè)置為3。但此分割結(jié)果不僅有正常細(xì)胞，還包含染色劑污染點(diǎn)，甚至有的淺染色細(xì)胞并沒有分割出來。第二步，算法使用Canny算子檢測細(xì)胞邊緣，經(jīng)過孔洞填充和腐蝕，重新確定了一組細(xì)胞核的位置坐標(biāo)。然后對屬于同一細(xì)胞核、被重復(fù)檢出的細(xì)胞位置進(jìn)行歸并，最終用細(xì)胞的特征參數(shù)排除染色劑污染點(diǎn)，提高了算法精度。

2 細(xì)胞圖像的K-means聚類分割

K-means算法是一種動(dòng)態(tài)聚類算法[12]，以歐式距離作為相似度測度，以誤差的平方和為聚類準(zhǔn)則函數(shù)，迭代計(jì)算使得準(zhǔn)則函數(shù)至收斂為止。算法目的是把n個(gè)樣本點(diǎn)分為k個(gè)簇，使簇內(nèi)具有較高的相似度。算法先隨機(jī)地選取k個(gè)對象作為初始聚類中心(聚類中心代表簇的平均值)，然后對剩余樣本根據(jù)到各個(gè)聚類中心的距離，將它們分配給最近的簇，再重新計(jì)算每個(gè)簇的平均值。此過程不斷迭代直至準(zhǔn)則函數(shù)J收斂。若進(jìn)行第m次迭代，需先更新第i個(gè)聚類中心zi

(1)

再計(jì)算準(zhǔn)則函數(shù)

(2)

其中，聚類中心數(shù)為k，Ni表示第i個(gè)簇的樣本數(shù)。本文設(shè)置聚類中心數(shù)為3。原圖像灰度化后，其像素被分成了3個(gè)簇(即細(xì)胞核、細(xì)胞質(zhì)、環(huán)境背景)，所以整張圖也被分割成了三值圖像。提取圖像中的細(xì)胞核，并用細(xì)胞核的質(zhì)心坐標(biāo)代表每個(gè)細(xì)胞的位置，最終得到位置集合S1(I)。此時(shí)，S1(I)可能包含染色劑污染點(diǎn)，也可能遺漏染色較淺導(dǎo)致未被分割的細(xì)胞核。

3 Canny算子邊緣檢測分割

Canny算子是基于最優(yōu)化算法的邊緣檢測算子[13]。實(shí)驗(yàn)證明，Canny算子在處理高斯白噪聲污染的圖像方面優(yōu)于其它傳統(tǒng)的邊緣算子。Canny 算子的實(shí)現(xiàn)主要包括四個(gè)部分：1)平滑圖像；2)計(jì)算梯度的幅值和方向；3)梯度方向上對梯度幅值做非極大值抑制；4)雙閾值方法檢測圖像邊緣。圖像中，每個(gè)像素的灰度的梯度值為

(3)

梯度方向?yàn)?/p>

H(x,y)=arctan (kx(x,y),ky(x,y))

(4)

其中，kx和ky分別為像素(x,y)的水平梯度和垂直梯度。

Canny算子作用后，圖像I生成了關(guān)于細(xì)胞邊緣的二值圖像。使用形態(tài)學(xué)進(jìn)行孔洞填充和腐蝕去噪，就得到了細(xì)胞核的質(zhì)心坐標(biāo)。這樣，圖像I產(chǎn)生了關(guān)于細(xì)胞位置的集合S2(I)。這里S2(I)包含了未被K-means檢出的細(xì)胞位置。而屬同一細(xì)胞核的質(zhì)心位置一定會(huì)被重復(fù)檢出，需要?dú)w并。所以算法需要考查上述二種分割得到位置的并集

S(I)=S1(I)∪S2(I)

(5)

4 結(jié)合Canny算子的必要性

(a)原圖；(b)K-means；(c)Canny 圖2 染色不足的情況舉例Fig.2 The case of inadequate dyeing

關(guān)于染色問題的一種情況是細(xì)胞核染色不足(圖2-a)。單獨(dú)使用K-means聚類處理圖像，形成了3個(gè)聚類中心，它們一般代表了圖像中細(xì)胞核、細(xì)胞質(zhì)、環(huán)境背景3種區(qū)域的灰度平均值。K-means聚類后，算法是根據(jù)分離出的細(xì)胞核來統(tǒng)計(jì)細(xì)胞數(shù)量，若染色不足(圖2-b)，該細(xì)胞是不會(huì)被統(tǒng)計(jì)的。但實(shí)踐中，即使染色不足導(dǎo)致細(xì)胞核灰度接近細(xì)胞質(zhì)，聚類方法失效，也還是能用Canny算子檢測出該細(xì)胞核的邊緣(圖2-c)，彌補(bǔ)K-means算法的缺陷。

(a)原圖；(b)K-means；(c)Canny 圖3 污染塊與正常細(xì)胞邊界模糊舉例Fig.3 The fuzzy boundary between the contaminant and normal cell

另一種情況是染色污染塊過于靠近正常細(xì)胞，兩者邊界模糊(圖3-a)。K-means聚類之后，污染塊和細(xì)胞核連成一體(圖3-b)，導(dǎo)致要識別的細(xì)胞核像素面積遠(yuǎn)超均值，被算法當(dāng)作大污染塊自動(dòng)排除了，此為統(tǒng)計(jì)錯(cuò)誤。但由于污染塊的顏色常不均勻，Canny算子作用后，污染塊不易形成閉合邊界(圖3-c)，孔洞填充后就不易形成實(shí)心點(diǎn)。這樣污染塊就和正常的細(xì)胞分開了，也彌補(bǔ)了K-means算法的缺陷。

如果單獨(dú)使用Canny算子分割圖像，一方面會(huì)遺漏沒有形成閉合邊界的細(xì)胞，即對這樣的細(xì)胞Canny算子是失效的。另一方面，也無法得到細(xì)胞核、細(xì)胞質(zhì)的像素面積、周長等量化特征為算法的后續(xù)處理做準(zhǔn)備。綜上所述，Canny算子是K-means聚類的必要補(bǔ)充。

5 染色劑污染點(diǎn)的排除

當(dāng)染色污染面積遠(yuǎn)大于細(xì)胞的平均值時(shí)，可以通過像素面積直接排除；當(dāng)污染面積與細(xì)胞大小相當(dāng)時(shí)，可以綜合細(xì)胞形狀特性參數(shù)來排除。第一種情況較為簡單，這里不再贅述，主要討論第二種情況。

圖4 像素計(jì)算范圍Fig.4 Pixel computing range

6 實(shí)驗(yàn)

下面用二例細(xì)胞顯微圖像測試所提出方法的有效性，并把結(jié)果與其它算法做了對比，顯示新算法的優(yōu)越性。

6.1 算例1

圖5 算例1的細(xì)胞圖像Fig.5 Original cell image of ex.1

圖5是一張經(jīng)過瑞士吉姆薩染色的原細(xì)胞圖像，其中有淺染色細(xì)胞(圖5藍(lán)虛線框)和若干染色劑污染點(diǎn)(圖5紅實(shí)線框)。

1)算法的有效性

第一次分割使用聚類中心數(shù)為3的K-means法得到三值圖像，如圖6所示。在圖6中提取白色的細(xì)胞核，再經(jīng)孔洞填充和腐蝕去噪，得到二值圖像(圖7)。圖7中發(fā)現(xiàn)，K-means分割出了大部分細(xì)胞核，但包含了不該包含的染色劑污染點(diǎn)(圖7紅實(shí)線框)，卻沒包含應(yīng)該包含的淺染色細(xì)胞(圖7白虛線框)。

然后使用Canny算子做二次分割獲取細(xì)胞邊緣特征，如圖8所示。經(jīng)孔洞填充(圖9)，再腐蝕去噪，由圖10計(jì)算出細(xì)胞核的質(zhì)心坐標(biāo)。分割效果令人滿意，因?yàn)閳D5中的藍(lán)虛線框所標(biāo)識的淺染色細(xì)胞，在圖9或圖10的對應(yīng)地方被Canny算子檢測出來了。

求圖7、圖10的位置的并集，得到圖11。其中紅圈表示K-means聚類得到的細(xì)胞位置，藍(lán)叉表示Canny算子得到的細(xì)胞位置。對屬于同一細(xì)胞核、被重復(fù)檢出的細(xì)胞質(zhì)心位置進(jìn)行歸并，并利用細(xì)胞核占比、細(xì)胞核面積等參數(shù)綜合篩選出真實(shí)的細(xì)胞。對于此例，設(shè)置細(xì)胞核面積閾值為22像素、核占比取值為0.25

圖6 K-means分割Fig.6 K-means segmentation

圖7 提取細(xì)胞核并腐蝕Fig.7 Nuclear extraction and corrosion

圖8 Canny算子檢測邊緣

圖9 孔洞填充Fig.9 Hole filling

圖10 腐蝕和去噪Fig.10 Corrosion and denoising

圖11 合并位置標(biāo)識Fig.11 Merge position identification

2)與其它方法比較

本文提出的方法與單獨(dú)使用K-means方法、Canny算子法、脈沖耦合神經(jīng)網(wǎng)絡(luò)(PCNN)分割和人工統(tǒng)計(jì)的方法對比，統(tǒng)計(jì)相對誤差如表1所示。K-means方法結(jié)果如圖7所示，Canny算子法標(biāo)識如圖9所示，PCNN神經(jīng)網(wǎng)絡(luò)標(biāo)識如圖13所示。從表1可以看出，本文提出的方法與人工統(tǒng)計(jì)相比，相對誤差最小(1.1%)，精度達(dá)到了98.9%。

圖12 合并標(biāo)識并排除染色劑污染點(diǎn)Fig.12 Eliminate contamination points

圖13 PCNN分割后再腐蝕Fig.13 PCNN and corrosion

6.2 算例2

用另一幅經(jīng)瑞士染色的細(xì)胞圖像來驗(yàn)證新方法的有效性，如圖14所示。其中有淺染色細(xì)胞(藍(lán)虛線框)和若干染色劑污染點(diǎn)(紅實(shí)線框)。

1)算法的有效性

經(jīng)中值濾波且使用聚類中心為3的K-means法分割，提取細(xì)胞核并去噪得到二值圖像，如圖15所示。與上例不同，本例還需用分水嶺法分割粘連細(xì)胞(紅框中)，得到圖16。

表1 不同算法的統(tǒng)計(jì)精度Table 1 Statistical accuracy of different algorithms

計(jì)算圖16中所有細(xì)胞核的質(zhì)心位置，結(jié)果標(biāo)注在圖17上。

圖14 算例2的細(xì)胞圖像Fig.14 Original cell image of ex.2

圖15 提取細(xì)胞核

圖16 分割粘連細(xì)胞Fig.16 Clustered cells separated

圖17 K-means法標(biāo)注的細(xì)胞位置Fig.17 Cell locations labeled by the K-means

以上看出，K-means法標(biāo)出了大部分細(xì)胞核的位置，但也存在著問題。17號點(diǎn)是染色污染點(diǎn)已被K-means排除，而 45和54號染色污染塊并沒有被K-means排除，這2點(diǎn)為錯(cuò)誤統(tǒng)計(jì)；另外， 6號點(diǎn)是細(xì)胞核(藍(lán)虛線框)，也沒能被K-means法識別。原因是K-means無法區(qū)分與6號點(diǎn)緊密相連的染色劑污染塊(見圖16)，導(dǎo)致提取的細(xì)胞核面積遠(yuǎn)超平均值被算法自動(dòng)排除了。限于篇幅，類似問題不一一敘述。

所以必須借助Canny算子做圖像的二次分割，輔助統(tǒng)計(jì)。Canny算子獲取細(xì)胞邊緣特征后，經(jīng)過孔洞填充(圖18)、分水嶺分割粘連細(xì)胞和腐蝕去噪(圖19)等步驟獲得細(xì)胞位置(在圖18中用藍(lán)×標(biāo)注)。將圖18結(jié)果與圖17對比，可以看出原6號位置的細(xì)胞在圖18中與染色劑污染塊分開了；且原54號位置的染色污染塊在圖18中也被Canny算子自動(dòng)排除了。

求圖17、18的位置的并集，即對屬于同一細(xì)胞核、被重復(fù)檢出的細(xì)胞質(zhì)心位置進(jìn)行歸并，并利用細(xì)胞核占比、核面積等參數(shù)綜合篩選細(xì)胞，最終得到圖20所示結(jié)果(統(tǒng)一用紅圈標(biāo)識位置)。從圖20看出，新方法排除了染色劑污染的干擾，提高了統(tǒng)計(jì)精度。

圖18 Canny算子檢測Fig.18 Canny operator detection

圖19 腐蝕和去噪Fig.19 Corrosion and denoising

圖20 合并標(biāo)識并排除染色劑污染點(diǎn)Fig.20 Eliminate contamination points

2)與其它方法比較

算例2與單獨(dú)使用K-means方法、Canny算子法、脈沖耦合神經(jīng)網(wǎng)絡(luò)(PCNN)分割和人工統(tǒng)計(jì)的方法對比，統(tǒng)計(jì)相對誤差如表2所示。可以看出，本文的方法與人工統(tǒng)計(jì)相比，相對誤差最小(2.3%)，精度達(dá)到了97.7%。

表2 不同算法的統(tǒng)計(jì)精度Table 2 Statistical accuracy of different algorithms

7 結(jié)論

本文提出了基于K-means聚類與Canny算子相結(jié)合方法用于細(xì)胞顯微圖像分割和統(tǒng)計(jì)。對于染色程度較淺、或有其它染色污染導(dǎo)致K-means聚類法未能正常識別細(xì)胞核的情況，本文采用Canny算子輔助分割以提高統(tǒng)計(jì)的準(zhǔn)確性。此方法不僅降低了細(xì)胞分析中對圖像染色質(zhì)量的要求，還有效地解決了現(xiàn)有細(xì)胞統(tǒng)計(jì)方法的誤差較大的問題。