李衡，龐楓騫，姜彬，劉志文

(1.北京理工大學(xué) 信息與電子學(xué)院，北京 100081; 2.北京理工大學(xué) 光電學(xué)院，北京 100081;3.國家留學(xué)基金管理委員會 信息資源部，北京 100081)

細(xì)胞是生物體最小的組成單元，其狀態(tài)能夠直接反映生物體的許多生理信息，因此細(xì)胞常被作為重要的研究對象，用于生物醫(yī)學(xué)領(lǐng)域的研究工作. 另一方面，細(xì)胞的形態(tài)受到細(xì)胞自身生理活動和所處環(huán)境的共同影響，能夠提供多種細(xì)胞生理狀態(tài)信息，是細(xì)胞篩選、識別和分類的重要依據(jù)[1-2].

近年來，隨著顯微成像技術(shù)的發(fā)展，細(xì)胞圖像序列被越來越多的用于細(xì)胞生理研究中，細(xì)胞形態(tài)的變化也逐漸成為研究者觀察的重要目標(biāo). 細(xì)胞的動態(tài)形變能夠更快速、準(zhǔn)確地反映細(xì)胞的狀態(tài)和活動，對檢測細(xì)胞狀態(tài)、預(yù)測細(xì)胞周期、研究細(xì)胞生理具有重要意義[3-4]. 卡耐基梅隆大學(xué)Murphy教授領(lǐng)導(dǎo)的團(tuán)隊(duì)[5]在利用形態(tài)特征描述細(xì)胞內(nèi)蛋白質(zhì)分布狀態(tài)的研究方面不斷有新的進(jìn)展，報道了亞細(xì)胞局部特征提取，使用三維圖像分析蛋白質(zhì)分布[6]，將蛋白質(zhì)分布分析應(yīng)用于癌癥早期檢測[7]等多項(xiàng)工作. 南開大學(xué)趙新教授帶領(lǐng)研究組對使用細(xì)胞形態(tài)檢測細(xì)胞活性的方法進(jìn)行了研究，提出一種基于形態(tài)參數(shù)度量細(xì)胞活性的方法[8]. 因此準(zhǔn)確高效地度量和分類細(xì)胞動態(tài)形變能夠?yàn)榧?xì)胞生物學(xué)研究提供強(qiáng)力幫助.

本文以植皮后小鼠淋巴細(xì)胞的形變作為研究對象，根據(jù)淋巴細(xì)胞形變的特性，設(shè)計(jì)出一種能夠有效處理細(xì)胞形變過程中空時信息的特征提取方法，并將其應(yīng)用于兩個淋巴細(xì)胞形變數(shù)據(jù)集上. 實(shí)驗(yàn)結(jié)果表明，所提出的方法能夠準(zhǔn)確度量和分析淋巴細(xì)胞的動態(tài)形變，相比于其他形變度量算法具有更高的準(zhǔn)確性和更好的魯棒性.

1 淋巴細(xì)胞形變數(shù)據(jù)

淋巴細(xì)胞的形變數(shù)據(jù)由合作單位，北京佑安醫(yī)院負(fù)責(zé)采集. 淋巴細(xì)胞采集自經(jīng)過植皮手術(shù)的小鼠(周齡6～8，體重20～22 g)外周血，使用相差顯微鏡(奧林帕斯BX51，放大倍率1 000×)進(jìn)行觀察，每次僅觀察一個淋巴細(xì)胞. 將目標(biāo)淋巴細(xì)胞置于顯微鏡視野中心后，并使用超顯微真彩圖像分析系統(tǒng)對影像進(jìn)行進(jìn)一步放大和拍攝. 以25幀/s的速度拍攝淋巴細(xì)胞20～30 s，保存為分辨率288×352的AVI格式視頻數(shù)據(jù). 隨后，對視頻數(shù)據(jù)進(jìn)行篩選，剔除拍攝不清晰，或淋巴細(xì)胞與其他細(xì)胞有重疊的數(shù)據(jù). 圖1展示了視頻數(shù)據(jù)的截圖，其中每行對應(yīng)一個視頻的圖像序列，每列圖片采樣時間間隔為5 s，完整的視頻數(shù)據(jù)示例可從http://isip.bit.edu.cn/kyxz/xzlw/77051.htm下載.

數(shù)據(jù)采集中，分別使用兩種小鼠實(shí)驗(yàn)方式，獲得兩個淋巴細(xì)胞形變數(shù)據(jù)集.

數(shù)據(jù)集Ⅰ：包含兩個小鼠植皮實(shí)驗(yàn)組，第一組使用20個健康雄性Balb/C小鼠進(jìn)行自體皮膚移植；第二組使用20個健康雄性Balb/C小鼠作為受體，20個健康雄性C57BL/6小鼠作為皮膚供體，進(jìn)行異體皮膚移植. 14 d后，斷尾采集小鼠外周血，觀察和拍攝淋巴細(xì)胞形變數(shù)據(jù). 將組別作為數(shù)據(jù)標(biāo)簽，每組各包含40個視頻數(shù)據(jù).

數(shù)據(jù)集Ⅱ：僅有一個小鼠植皮實(shí)驗(yàn)組，使用25個健康雄性Balb/C小鼠作為受體，25個健康雄性C57BL/6小鼠作為皮膚供體，進(jìn)行異體皮膚移植. 7 d后，斷尾采集小鼠外周血，觀察和拍攝淋巴細(xì)胞形變數(shù)據(jù). 由3名經(jīng)過訓(xùn)練的觀察者對淋巴細(xì)胞動態(tài)形變進(jìn)行判別，將視頻數(shù)據(jù)人工分為4組(正常組，輕度活躍組，中度活躍組和高度活躍組)，最終通過投票方式確定數(shù)據(jù)標(biāo)簽，每組各包含25個視頻數(shù)據(jù).

2 形變局部時序特征提取

通過對細(xì)胞視頻數(shù)據(jù)的詳細(xì)觀察，可以發(fā)現(xiàn)細(xì)胞的形態(tài)變化不僅在時間維度上有連續(xù)性，還在空間的各個局部區(qū)域間表現(xiàn)出不一致性. 為了更直觀地展示細(xì)胞形變的特性，圖2中給出了數(shù)據(jù)集I中兩類細(xì)胞動態(tài)形變的對比圖. 圖中已提取出細(xì)胞圖像序列的輪廓，并使用細(xì)胞輪廓序列繪制出細(xì)胞的動態(tài)形變. 可以看出，對于正常形變的細(xì)胞，其輪廓相對穩(wěn)定，整體變化不明顯，只在少數(shù)區(qū)域能看到輪廓的演化. 與之相對，異常形變細(xì)胞的輪廓變化顯著，在形變劇烈的區(qū)域還可以看出較明顯的連續(xù)性；同時整個細(xì)胞各部分形變程度并不統(tǒng)一，圖中左側(cè)區(qū)

域的形變程度遠(yuǎn)高于右側(cè). 根據(jù)觀察結(jié)果，設(shè)計(jì)出一種能夠有效處理淋巴細(xì)胞形變特性的特征提取算法.

2.1 基于徑向距離的細(xì)胞動態(tài)形變描述

細(xì)胞的輪廓變化能夠很好地反映其形變過程，因此首先利用細(xì)胞輪廓序列描述細(xì)胞的整個形變過程. 考慮到這里所觀察的淋巴細(xì)胞主要呈圓形和橢圓形，具體使用圖像特征中的徑向距離方法，表示細(xì)胞輪廓序列. 將在角度k處，時間點(diǎn)t的輪廓點(diǎn)的徑向距離定義為

k∈[0,360)，t=1,2,…,T.

(1)

式中：細(xì)胞圖像序列第一幀的細(xì)胞質(zhì)心(xc,yc)被作為提取徑向距離的參考點(diǎn);(xk,t,yk,t)表示角度k處，時間點(diǎn)T的輪廓點(diǎn)位置;T為細(xì)胞圖像序列的長度. 圖3(a)3(b)中展示了使用徑向距離表示的細(xì)胞形變輪廓序列.

隨后，從時間維度分析細(xì)胞輪廓點(diǎn)的形變信息，將式(1)中的徑向距離函數(shù)表示成細(xì)胞輪廓上角度k處隨時間t的變化函數(shù)Ck(t). 如圖3(c)3(e)所示，以角度為120°的輪廓點(diǎn)為例，將細(xì)胞輪廓的徑向距離函數(shù)表示為其隨時間t的變化，記為C120(t)，其曲線如圖3(c)3(e)中的灰色曲線.

從圖中可以看出，細(xì)胞輪廓點(diǎn)的徑向距離函數(shù)在其隨時間t變化的過程中趨勢并不穩(wěn)定，灰色曲線上有許多抖動. 對圖中的灰色曲線進(jìn)行比較，在相對較短的時間段內(nèi)(大約20幀以內(nèi))，受到曲線抖動的影響，人眼無法區(qū)分圖中灰色曲線的差別. 幸運(yùn)的是，當(dāng)從一個較長的時間段觀察灰色曲線，圖3(e)中的曲線有比較明顯的變化趨勢.

根據(jù)這一觀察結(jié)論，提出利用曲線平滑的方法增強(qiáng)該輪廓點(diǎn)變化的趨勢信息. 實(shí)際中，使用7階多項(xiàng)式平滑的方式，得到平滑后的曲線fk(t)，如圖3(c)3(e)中黑色曲線所示. 通過黑色曲線，能夠更清楚地觀察和區(qū)分兩個細(xì)胞輪廓點(diǎn)的不同變化過程.

2.2 局部時序特征

根據(jù)細(xì)胞輪廓點(diǎn)隨時間的變化過程，可以看出細(xì)胞的局部變化特性能夠反應(yīng)出細(xì)胞動態(tài)形變的重要信息. 因此，通過提取細(xì)胞變化的局部時序特征，能夠更加有效地度量細(xì)胞的動態(tài)形變.

以平滑后的細(xì)胞輪廓點(diǎn)隨時間的形變曲線為基礎(chǔ)，提取細(xì)胞輪廓點(diǎn)變化的時序特征. 根據(jù)觀察，細(xì)胞輪廓點(diǎn)的變化趨勢是整個變化過程中最顯著的特征. 因此，對平滑曲線fk(t)求導(dǎo)，得到其導(dǎo)數(shù)f′k(t)，如圖3(d)3(f)中灰色曲線所示. 鑒于細(xì)胞輪廓點(diǎn)的變化僅有遠(yuǎn)離中心和朝向中心點(diǎn)兩個方向，使用細(xì)胞輪廓點(diǎn)隨時間運(yùn)動的方向?qū)π巫冓厔葸M(jìn)行區(qū)分. 對應(yīng)到曲線f′k(t)，作為平滑曲線fk(t)的導(dǎo)數(shù)，其正值對應(yīng)遠(yuǎn)離中心的運(yùn)動，負(fù)值對應(yīng)朝向中心的運(yùn)動. 所以，可以使用f′k(t)取值的正負(fù)號區(qū)分輪廓點(diǎn)的運(yùn)動趨勢. 在同一變化趨勢下，計(jì)算輪廓點(diǎn)運(yùn)動的累積量Gk描述對應(yīng)形變的時序特征，如圖3(d)3(f)中的青色區(qū)域. 在實(shí)際計(jì)算時，以f′k(t)取值的符號變化點(diǎn)作為邊界，對取絕對值后的f′k(t)進(jìn)行積分，計(jì)算Gk，其定義為

(2)

式中：τj表示f′k(t)取值符號變化的點(diǎn)；T為細(xì)胞圖像序列的長度；j和J分別表示符合變化點(diǎn)的序數(shù)和總數(shù).

最后，計(jì)算Gk的最大值作為角度k處細(xì)胞輪廓點(diǎn)的局部時序變化量，并對整個數(shù)據(jù)集中的細(xì)胞各輪廓點(diǎn)的變化量進(jìn)行統(tǒng)計(jì)，繪制20個平均間隔的直方圖. 將某個細(xì)胞圖像序列在該直方圖各區(qū)間的分布作為其局部時序形變特征向量.

3 實(shí)驗(yàn)結(jié)果及分析

使用淋巴細(xì)胞形變數(shù)據(jù)庫驗(yàn)證所提算法的效果，并與若干個使用圖像特征的形變度量算法進(jìn)行比較，分析所提算法的效果. 實(shí)驗(yàn)中，使用高斯核函數(shù)的支持向量機(jī)(SVM)作為分類器訓(xùn)練分類模型. 整個分類過程使用5折交叉驗(yàn)證，獨(dú)立隨機(jī)重復(fù)100次，最終使用精確率(真陽性/(真陽性+假陽性))，召回率(真陽性/(真陽性+假陰性))和F值((2×精確率×召回率)/(精確率+召回率))作為評價分類效果的標(biāo)準(zhǔn).

實(shí)驗(yàn)數(shù)據(jù)中的淋巴細(xì)胞預(yù)先使用主動輪廓模型進(jìn)行了分割， 進(jìn)行了尺度歸一化和序列對齊的操作. 這里分別使用徑向距離、傅里葉描述子[9]、通用傅里葉描述子(generic fourier descriptor，GFD)[10]、Zernike矩[11]和閾值鄰域統(tǒng)計(jì)(threshold adjacency statistics，TAS)[12]，5種圖像特征為基礎(chǔ)的算法分析細(xì)胞形變，并進(jìn)行比較. 此外，前期工作提出的匹配游程統(tǒng)計(jì)方法的結(jié)果也一起給出. 具體計(jì)算時，徑向距離和傅里葉描述子都取180維，通用傅里葉描述子的角度和半徑分辨率為60和180， Zerniek矩共30階. 以1 s(25幀)為間隔對細(xì)胞視頻采樣后，對樣本圖像序列依次計(jì)算特征值，隨后使用歐氏距離量化相鄰樣本圖像間的差異. 最后，統(tǒng)計(jì)歐氏距離的值，并對全部數(shù)據(jù)的統(tǒng)計(jì)結(jié)果計(jì)算有20個平均間隔的直方圖，將其作為描述細(xì)胞形變的特征矢量，使用SVM分類.

3.1 數(shù)據(jù)集Ⅰ

首先在數(shù)據(jù)集Ⅰ上對現(xiàn)有的方法進(jìn)行驗(yàn)證. 數(shù)據(jù)集Ⅰ包含兩個分類，各有40個細(xì)胞視頻數(shù)據(jù)，5折交叉驗(yàn)證中，每類隨機(jī)選取32個數(shù)據(jù)為訓(xùn)練集，8個為測試集. 分類結(jié)果在訓(xùn)練集和測試集上的平均精確率(TrPr，TePr)，召回率(TrRe，TeRe)和F值(TrF，TeF)均已在表1中給出.

表1 數(shù)據(jù)集Ⅰ細(xì)胞形變分類結(jié)果Tab.1 Classification results of cell deformation in Dataset Ⅰ %

由于所研究的淋巴細(xì)胞一般為圓形或橢圓形，使用徑向距離能夠很好地提取細(xì)胞形變特征，因此其能夠表現(xiàn)出良好的效果. 傅里葉描述子是通過對目標(biāo)的輪廓曲線做一維傅里葉變化得到的，因此相比于其他特征提取方法，其對噪聲更加敏感，因此傅里葉描述子表現(xiàn)并不突出. 通用傅里葉描述子和Zernike矩都從目標(biāo)圖像中提取高維特征矢量描述形狀，能夠分辨不同形狀間的細(xì)小差異，是處理復(fù)雜形狀的有力工具，在細(xì)胞形態(tài)研究中有廣泛應(yīng)用[5, 13-14]. 在研究細(xì)胞形變時，使用通用傅里葉描述子和Zernike矩得到的特征矢量可以有效描述細(xì)胞形態(tài)的細(xì)節(jié)，以此為基礎(chǔ)計(jì)算出的特征空間距離能夠提取細(xì)胞形態(tài)間的細(xì)小差異，準(zhǔn)確描述形變特性. 閾值鄰域統(tǒng)計(jì)是用于描述細(xì)胞紋理的特征，其主要用于研究細(xì)胞區(qū)域內(nèi)的灰度分布變化. 雖然在細(xì)胞形變的過程中，細(xì)胞區(qū)域內(nèi)的紋理也會隨之變化，但是這種變化不是細(xì)胞形變的主要信息. 故而，基于紋理特征分析細(xì)胞形變無法發(fā)揮良好效果.

使用所提出的算法，細(xì)胞形變過程中的時序信息被最大程度地保存下來，有效提高了對形變的分類結(jié)果. 不僅如此，局部時序特征充分考慮的細(xì)胞形變中各個局部區(qū)域變化不一致的情況，分別度量細(xì)胞各輪廓點(diǎn)的時序形變，最后對細(xì)胞整體形變做出綜合評價. 所以，局部時序特征非常適合用于淋巴細(xì)胞動態(tài)形變的研究.

3.2 數(shù)據(jù)集Ⅱ

數(shù)據(jù)集Ⅱ包含4個分類，各有25個細(xì)胞視頻數(shù)據(jù)，使用5折交叉驗(yàn)證時，從各類別中隨機(jī)選取20個數(shù)據(jù)為訓(xùn)練集，5個為測試集. 表2中給出了各算法在數(shù)據(jù)集Ⅱ上的分類結(jié)果.

表2 數(shù)據(jù)集Ⅱ細(xì)胞形變分類結(jié)果Tab.1 Classification results of cell deformation in Dataset Ⅱ %

相比于數(shù)據(jù)集Ⅰ，各算法在數(shù)據(jù)集Ⅱ上的效果均明顯降低. 這首先是由于數(shù)據(jù)集Ⅱ上的多分類問題明顯比數(shù)據(jù)集Ⅰ上的二分類問題復(fù)雜，會影響算法的表現(xiàn). 其次，數(shù)據(jù)集Ⅰ中的兩個類別分別是從自體皮膚移植和異體皮膚移植小鼠的淋巴細(xì)胞采集而來，兩個類別間細(xì)胞形變差異顯著，比較易于分辨；而數(shù)據(jù)集Ⅱ中各類別是根據(jù)人工觀察標(biāo)記而來，細(xì)胞本身所處生理?xiàng)l件上沒有差異，形變間區(qū)分度相對較低. 而且受到人工標(biāo)記影響，數(shù)據(jù)集Ⅱ的分類標(biāo)簽中包含人為偏差，更提高了分類的難度.

在分類效果上，對比算法的效果下降比所提出的算法更加顯著，這是由于對數(shù)據(jù)集Ⅱ中的細(xì)胞形變進(jìn)行分類更加需要充分提取形變信息. 直接從相鄰幀的形態(tài)差異提取形變，對比算法對圖像序列進(jìn)行采樣，再從樣本幀間提取形變信息的框架本身就丟失了許多未被采樣到的形變信息. 然而，由于顯微視頻數(shù)據(jù)幀間隔時間短，對全部幀提取形態(tài)，信噪比低，難以從噪聲中分辨有用信息.

為了更有效地展示所提算法的效果，將數(shù)據(jù)集Ⅱ上每個視頻數(shù)據(jù)對應(yīng)的局部時序特征直方圖在圖4中進(jìn)行展示. 圖中可以看出，細(xì)胞形變正常組的特征值整體較小，且分布比較集中. 而隨著細(xì)胞形變活躍度的提高，特征值不斷變大，且分布越來越分散. 此外，相鄰活躍度的組間，特征分布的區(qū)分度較低，分類邊界不清晰. 尤其輕度活躍組和中度活躍組的特征分布相似度非常高.

對各組間的局部時序特征計(jì)算巴氏距離(局部時序特征)，其定義為

(3)

式中：H為數(shù)值分布；γ為數(shù)值區(qū)間. 圖5中展示了各組間細(xì)胞形變局部時序特征的平均巴氏距離.

從圖5中可以看出，輕微活躍組和中度活躍組間的巴氏距離較大，說明兩組中的特征分布情況相似. 高度活躍組的組內(nèi)巴氏距離相比于其他組較小，表明該組特征分布的組內(nèi)相似度低. 由此看出，數(shù)據(jù)集Ⅱ的細(xì)胞動態(tài)形變多分類問題相對復(fù)雜，僅使用圖像特征度量全局形變難以有效解決此問題. 通過提取形變過程中的時序信息和空間信息，所提特征算法能更好地度量和分類細(xì)胞動態(tài)形變.

4 結(jié) 論

針對顯微視頻中細(xì)胞的動態(tài)形變進(jìn)行研究，根據(jù)對細(xì)胞形變的觀察，提出一種針對細(xì)胞形變的局部時序特征，能夠有效提取細(xì)胞形變過程中時序信息和空間信息. 使用兩個淋巴細(xì)胞形變視頻數(shù)據(jù)集，對該特征算法的效果進(jìn)行了驗(yàn)證. 實(shí)驗(yàn)結(jié)果表明，所提特征算法在兩個數(shù)據(jù)集上均有突出表現(xiàn)，能夠準(zhǔn)確度量和分類細(xì)胞動態(tài)形變，相比于基于圖像特征的形變分析算法具有更大優(yōu)勢.