基于傅里葉變換的細(xì)胞牽引力測量研究

溫 娜，張 帆,2

（1.北京信息科技大學(xué) 光電測試技術(shù)及儀器教育部重點(diǎn)實(shí)驗(yàn)室，北京 100192；2.北京信息科技大學(xué)先進(jìn)光電子器件與系統(tǒng)創(chuàng)新引智基地，北京 100192）

引言

生命從誕生開始就生活在力學(xué)環(huán)境中并與之相適應(yīng)。細(xì)胞與細(xì)胞外微環(huán)境的力學(xué)相互作用對于細(xì)胞的遷移、增殖和分化都起著至關(guān)重要的作用[1-4]。但由于細(xì)胞內(nèi)在的生物復(fù)雜性及微納米尺度的制約，傳統(tǒng)的宏觀實(shí)驗(yàn)力學(xué)方法與技術(shù)往往難以直接使用。經(jīng)微加工技術(shù)制備的垂直排列的聚二甲基硅氧烷（ploydimethylsiloxane，PDMS），在頂端表面修飾后可支持細(xì)胞粘附，細(xì)胞粘附后產(chǎn)生的牽引力會造成微柱彎曲，將測得的撓度變化輸入簡單的材料力學(xué)公式即可獲得細(xì)胞的牽引力分布。采用這種致密、垂直、離散微柱矩陣結(jié)構(gòu)替代傳統(tǒng)測量的連續(xù)介質(zhì)，通過對微柱形變的顯微圖像處理，細(xì)胞牽引力可以被直接定量測量，精度可以達(dá)到數(shù)十nN·μm?1量級[5-6]。

微柱形變的顯微圖像處理，作為細(xì)胞牽引力測量的最后一環(huán)，直接決定著細(xì)胞牽引力最終測量結(jié)果的精度和準(zhǔn)確性。而在微柱的顯微圖像處理中，微柱頂端的位移識別及提取精度又是影響數(shù)字圖像處理的核心因素之一。

目前在數(shù)字圖像處理中常用的微柱頂端位移識別的方法有質(zhì)心法（centroid method）和霍夫變換法（Hough transform）等。質(zhì)心法利用微柱在高倍顯微成像下成像模型圓對稱性質(zhì)，利用Matlab中的計(jì)算區(qū)域描繪子函數(shù)regionprops 計(jì)算質(zhì)心的位置從而獲得微柱位置。具體操作是將微柱的高倍顯微圖轉(zhuǎn)換為灰度圖，然后將圖像進(jìn)行中值濾波、二值化等處理，最后通過計(jì)算得到微柱端面質(zhì)心[7-12]?；舴蜃儞Q是數(shù)字圖像處理中常用的一種特征提取技術(shù)。它將一個(gè)空間中具有相同形狀的曲線或直線映射到另一個(gè)坐標(biāo)空間的一個(gè)點(diǎn)上形成峰值，從而把檢測任意形狀的問題轉(zhuǎn)化為統(tǒng)計(jì)峰值問題。利用霍夫變換進(jìn)行微柱頂端的位移識別，實(shí)質(zhì)上是用霍夫變換提取分布于目標(biāo)圓周上的參數(shù)及點(diǎn)的特征值來檢測圓，對圖像上的每個(gè)點(diǎn)定義一個(gè)參數(shù)空間的映射，通過在參數(shù)空間找尋特征值（峰值或最大值）得到為與圖像空間中的目標(biāo)形狀。具體操作是將高倍顯微下的微柱圖轉(zhuǎn)化為灰度圖，進(jìn)行雙邊濾波后通過霍夫變換檢測圓，得到圓的圓心和半徑[13-15]。在用霍夫變換檢測圓時(shí)，需要一個(gè)三維數(shù)組用來進(jìn)行累加運(yùn)算，因此當(dāng)所要檢測的圓太多時(shí)，計(jì)算量和內(nèi)存占用量將會急劇增大，計(jì)算效率將會大大降低。而利用regionprops 函數(shù)中的centroid 參數(shù)提取質(zhì)心時(shí)，是通過一階中心矩陣計(jì)算出的質(zhì)心，無論是在計(jì)算量上，還是內(nèi)存占用量上，都遠(yuǎn)遠(yuǎn)小于霍夫變換檢測圓時(shí)的計(jì)算量和內(nèi)存占用量，計(jì)算效率遠(yuǎn)大于霍夫變換的計(jì)算效率[16-19]。但是用質(zhì)心法處理高倍顯微鏡下的PDMS 微柱圖像時(shí)，可能會由于微柱陣列自身基底不平或者照明不均導(dǎo)致處理圖像時(shí)很難選取到合適的閾值，進(jìn)而影響后續(xù)的一系列操作；而且質(zhì)心法處理圖像需要有原始位置作為標(biāo)準(zhǔn)，拍攝原始位置或者是算法的處理都較為繁瑣；此外，質(zhì)心法二值化的過程致使其輪廓跳動較大，會給微柱端面質(zhì)心定位帶來誤差。因此質(zhì)心法雖相較于霍夫變換法運(yùn)算速度較快，但受這種方法自身的局限性，仍存在著照明不均時(shí)閾值的合理選取困難，需原始位置作為標(biāo)準(zhǔn)，二值化的過程使得輪廓跳動大這些問題。

為了能更直觀、更快速地得到細(xì)胞牽引力的分布情況，克服質(zhì)心法處理過程中存在的一些問題，本文提出了一種基于數(shù)字傅里葉變換的細(xì)胞牽引力測量方法。該方法通過預(yù)設(shè)變形二維點(diǎn)陣圖像仿真實(shí)驗(yàn)，驗(yàn)證了測量方法的可行性，探究了所選濾波器種類、濾波窗口大小、偏移點(diǎn)大小對實(shí)驗(yàn)結(jié)果的影響。測量得到了微柱高倍顯微圖所受細(xì)胞力的分布情況并進(jìn)一步證明了該方法與質(zhì)心法、霍夫變化法相比在運(yùn)算速度和運(yùn)算結(jié)果呈現(xiàn)方面的優(yōu)越性。

1 測量原理

當(dāng)培養(yǎng)在微柱上的細(xì)胞的牽引力帶動微柱偏轉(zhuǎn)，從而使偏轉(zhuǎn)微柱部分的周期性微結(jié)構(gòu)發(fā)生改變時(shí)，可通過二維傅里葉變換解調(diào)出變形區(qū)域內(nèi)的各單元變形量。以二維光柵為例，假設(shè)一個(gè)未變形的理想正弦光柵，其柵線與X方向、Y方向均不平行。將該未變形的光柵用作參考光柵，其像面光強(qiáng)分布可以表示為

式中：fx為未變形的參考光柵在X方向上的頻率；fy為未變形的參考光柵在Y方向上的頻率； φ0為原始光柵像的初始相，包含變形的光柵，由于變形引起了相位調(diào)制 φ (x,y)；r(x,y)為待測物體表面非均勻的反射率，將變形柵的像光強(qiáng)分布展開為傅里葉級數(shù)，得到：

根據(jù)二維DFT 的可分離性改寫為

式中：

若未變形柵線與一軸線平行，如與Y軸平行，即fy=0，此時(shí)未變形柵線的光強(qiáng)分布的傅里葉級數(shù)可寫為

未變形柵線與Y軸平行，此時(shí)變形柵線的光強(qiáng)分布傅里葉級數(shù)可寫為

式中：r(x,y)為待測物體表面非均勻的反射率；φ(x,y)為由變形引起的相位調(diào)制。

（5）式、（6）式僅是未變形柵線與Y軸平行時(shí)的特例，對于其他更為一般的情況，對（1）式和（2）式兩式作二維傅里葉變換，并用濾波器將低頻直流號和高頻雜散信號濾掉，再將基頻部分進(jìn)行二維傅里葉逆變換，得到變形后柵線的光強(qiáng)分布：

對比未變形柵線的光強(qiáng)分布：

可知，當(dāng)光柵的局部發(fā)生變形時(shí)，變形處的周期性被破壞，此時(shí)變形處的頻率遠(yuǎn)小于其他未變形處的頻率，相當(dāng)于（7）式中的 Δfx、 Δfy都是負(fù)值。因此，對應(yīng)到像面復(fù)振幅分布中，發(fā)生變形的光柵其形變部位的峰值也會小于未發(fā)生形變部位的峰值，且變形程度越大對應(yīng)的復(fù)振幅分布中的振幅也越小。

光強(qiáng)到細(xì)胞牽引力的計(jì)算，依賴于光強(qiáng)到微柱偏移量與微柱偏移量到細(xì)胞力兩方面的理論模型。光強(qiáng)g和微柱位移 δ的對應(yīng)關(guān)系，可由Matlab作對微柱作二維光場仿真計(jì)算得出數(shù)學(xué)模型：

其中A、B兩值根據(jù)微柱直徑、排布方式以及間距來確定，本文中對直徑3 μm、間距9 μm、90°排布方式的微柱代入計(jì)算，得到A=50.793 4，B=?20.656 9，具體分析由仿真實(shí)驗(yàn)得出。

細(xì)胞牽引力F可由微柱位移量 δ得到：

式中：E為微柱的剛度，由COMSOL Multiphysics對微柱建立的數(shù)學(xué)模型得到。

2 實(shí)驗(yàn)方法

2.1 預(yù)設(shè)變形二維點(diǎn)陣仿真

首先利用Photoshop，參照實(shí)際中常見的微柱尺寸（直徑為3 μm，相鄰間距為9 μm），等比例繪制30×30 的二維點(diǎn)陣圖以及二維點(diǎn)陣中某點(diǎn)發(fā)生位移的點(diǎn)陣圖，如圖1（a）為未變形的二維點(diǎn)陣圖，圖1（b）、圖1（c）、圖1（d）為偏移點(diǎn)分別在中心、左上角、右下角的二維點(diǎn)陣圖。通過算法將點(diǎn)陣圖像填充為原圖像陣列大小的兩倍來進(jìn)行0 填充處理。再通過Matlab 軟件對預(yù)設(shè)變形的二維點(diǎn)陣圖形分別進(jìn)行二維快速傅里葉變換、選取特定的高頻點(diǎn)進(jìn)行濾波、逆傅里葉變換等處理得到其對應(yīng)的幅值分布圖，以探究二維點(diǎn)陣圖中單點(diǎn)發(fā)生位置偏移以及一些列點(diǎn)均發(fā)生偏移時(shí)得到的幅值圖中的結(jié)果。

圖1 預(yù)設(shè)偏移的二維點(diǎn)陣圖Fig.1 Two-dimensional lattice diagrams of default offset

預(yù)設(shè)變形的二維點(diǎn)陣圖（以圖1（c）為例）在經(jīng)二維快速傅里葉變換后，得到的頻譜圖及頻譜圖的局部放大圖如圖2 所示。選取其中光強(qiáng)較大，并且不包含低頻直流成分的二維點(diǎn)陣圖的一級衍射斑點(diǎn)分別進(jìn)行標(biāo)號，分別得到這8 個(gè)衍射斑點(diǎn)作傅里葉逆變換的幅值分布圖以探究單頻點(diǎn)的選擇對仿真實(shí)驗(yàn)結(jié)果的影響。在此基礎(chǔ)上將發(fā)生不同方向偏移的點(diǎn)置于同一圖中進(jìn)行相同的操作以直觀比較選取不同頻點(diǎn)對得到幅值圖的影響，以探究在偏移點(diǎn)的數(shù)量或其他信息未知時(shí)對頻點(diǎn)選取的最佳方案。

圖2 頻譜圖及其局部放大圖Fig.2 Spectrogram and its partial magnification diagram

如圖3，A、B、C、D分別為只發(fā)生X方向偏移的點(diǎn)、只發(fā)生Y方向偏移的點(diǎn)和發(fā)生右下、左上方向偏移的點(diǎn)，對這樣同時(shí)包含典型偏移方向的偏移點(diǎn)的圖像進(jìn)行傅里葉變換，選取不同頻點(diǎn)以相同的濾波窗口大小進(jìn)行濾波，探究頻點(diǎn)選擇對幅值圖的影響。除此之外，濾波器窗口大小的設(shè)置也直接影響仿真實(shí)驗(yàn)結(jié)果，因此以圖1（c）為例針對同一個(gè)頻點(diǎn)，只改變?yōu)V波器窗口的大小，探究不同濾波器及濾波器濾器窗口大小對得到幅值分布圖結(jié)果的影響。

圖3 不同方向偏移的點(diǎn)陣圖Fig.3 Lattice diagram with offset in different directions

由于實(shí)際中細(xì)胞黏附在微柱上，通過細(xì)胞牽引力帶動微柱偏轉(zhuǎn)，可能出現(xiàn)某個(gè)方向上細(xì)胞牽引力太大而使微柱產(chǎn)生較大的偏轉(zhuǎn)，此時(shí)的微柱的俯視圖中看到的可能只是該微柱頂端的投影，這樣得到的微柱直徑相比較于其他未偏轉(zhuǎn)的微柱會發(fā)生明顯改變。因此在前面仿真的基礎(chǔ)上，進(jìn)一步探究不同大小的偏移點(diǎn)對仿真實(shí)驗(yàn)結(jié)果的影響。如圖4 為預(yù)設(shè)某點(diǎn)發(fā)生偏移的二維點(diǎn)陣點(diǎn)陣圖的局部放大圖，紅框內(nèi)為大小可變的偏移點(diǎn)，其圓心坐標(biāo)為固定值。除偏移點(diǎn)外其余點(diǎn)的直徑都均為15 個(gè)像素。對不同大小的偏移點(diǎn)經(jīng)過傅里葉變換，選取同一高頻點(diǎn)以同一窗口大小對其進(jìn)行濾波，再逆傅里葉變換得到幅值分布圖以探究偏移點(diǎn)大小對仿真實(shí)驗(yàn)結(jié)果的影響。

圖4 偏移點(diǎn)陣局部放大圖Fig.4 Partial magnification diagram of offset lattice

2.2 細(xì)胞-微柱顯微照片數(shù)字圖像處理

在預(yù)設(shè)變形二維點(diǎn)陣圖的基礎(chǔ)上，用真實(shí)顯微鏡下拍到的高倍顯微圖替代前面仿真中用Photoshop繪制的二維點(diǎn)陣圖，對其進(jìn)行數(shù)字傅里葉變換、濾波、逆傅里葉變換等操作后，得到幅值分布圖，以進(jìn)一步探究利用數(shù)字傅里葉變換測量細(xì)胞牽引力的可行性。

具體的測量過程分為兩大部分，首先是微柱的制備、LX-2 肝星狀細(xì)胞的接種及后續(xù)培養(yǎng)，一段時(shí)間后在選取生長狀態(tài)良好且細(xì)胞對微柱產(chǎn)生了較明顯拉伸的細(xì)胞區(qū)域進(jìn)行拍照，得到60X 物鏡下的高倍顯顯微鏡微圖，如圖5 所示[20]。然后將高倍顯微圖轉(zhuǎn)化為灰度圖并填充為2 048×2 048 像素大小的圖像，再對其進(jìn)行二維快速傅里葉變換、濾波、傅里葉逆變換等處理得到微柱陣列的幅值分布圖，在濾波時(shí)選取高頻點(diǎn)1 進(jìn)行濾波，具體分析步驟如圖6 所示。

圖5 細(xì)胞高倍顯微圖Fig.5 High-power micrograph of cell

圖6 程序分析步驟圖Fig.6 Flow chart of process analysis

高倍顯微圖經(jīng)數(shù)字傅里葉變換后得到的是幅值分布圖，要想測得細(xì)胞力必須首先得到高倍顯微圖中微柱的偏移量，再代入力學(xué)模型得出細(xì)胞力。

為了將經(jīng)圖6 的程序分析得到的幅值圖轉(zhuǎn)換為微柱偏移量，需對偏移的值和對應(yīng)的幅值進(jìn)行標(biāo)定。具體標(biāo)定過程為：參照高倍顯微圖灰度圖中未發(fā)生偏移的微柱大小及與其他未偏移微柱之間的間距，構(gòu)建與高倍顯微圖經(jīng)填充后的灰度圖大小完全相同的灰度圖，以標(biāo)準(zhǔn)微柱直徑及間隔在圖上構(gòu)建17×17 個(gè)“微柱”陣列，除最中心一個(gè)“微柱”發(fā)生偏移外，其余皆呈周期性排布，如圖7 所示。

圖7 作為標(biāo)定的圖及其局部放大圖Fig.7 Diagram as calibration and its partial magnification diagram

設(shè)置中心偏移“微柱”的偏移量為從0 到60 像素，以10 像素為步長得到一系列的灰度圖。將這些灰度圖做為待測圖，選擇與細(xì)胞-微柱顯微圖的灰度圖進(jìn)行運(yùn)算時(shí)相同的參數(shù)重復(fù)圖6 所示程序運(yùn)算，得到一系列幅值圖，測量得到不同位移幅度下的幅值分布的系統(tǒng)響應(yīng)，如圖8 所示。標(biāo)定經(jīng)數(shù)字傅里葉變換、濾波、傅里葉逆變換處理后幅值分布和局部偏移移量的關(guān)系。

圖8 偏移量及幅值的對應(yīng)關(guān)系圖Fig.8 Corresponding diagram of offset and amplitude

由圖8，微柱偏移量在0～60 像素（對應(yīng)實(shí)際距離為3 μm）范圍內(nèi)，經(jīng)OriginPro 程序可分析擬合出偏移量與幅值二者存在的線性關(guān)系，即通過標(biāo)定可以從幅值分布圖解算出微柱偏移的像素值，再由細(xì)胞的高倍顯微圖像素值換算微柱實(shí)際偏移量的分布。依據(jù)前文提到的本實(shí)驗(yàn)中微柱陣列的固化參數(shù)（固化溫度65 ℃，固化時(shí)間40 h），查閱相關(guān)文獻(xiàn)[21]得到此時(shí)PDMS 的楊氏模量約為2.66 MPa。在已有的相關(guān)文獻(xiàn)基礎(chǔ)上，調(diào)整楊氏模量和最大施加的力等參數(shù)，用多物理場仿真軟件COMSOL Multiphysics 對微柱頂端受力發(fā)生形變進(jìn)行仿真，得到此實(shí)驗(yàn)中微柱的力學(xué)模型。圖9 為在微柱頂端施加150 nN 力時(shí)微柱頂端的偏移情況，此時(shí)微柱底端保持不變，只在微柱頂端發(fā)生偏移，且此時(shí)最大偏移量為3.52 μm。通過仿真可知，在微柱頂端施加0～150 nN 的力時(shí)，微柱頂端偏移量與所受的力有很好的線性關(guān)系，且此時(shí)微柱剛度約為42.66 nN·μm?1，細(xì)胞牽引力分布情況最終由微柱的偏移量乘微柱剛度得到。

圖9 COMSOL Multiphysics 仿真結(jié)果Fig.9 Simulation results of COMSOL Multiphysics

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

3.1 預(yù)設(shè)變形二維點(diǎn)陣

3.1.1 預(yù)設(shè)變形點(diǎn)的位置

分別將圖1（a）、圖1（b）、圖1（c）、圖1（d）經(jīng)0 填充后進(jìn)行快速傅里葉變換，單頻點(diǎn)濾波，再逆傅里葉變換，得到的幅值分布圖如圖10 所示。

圖10 不同位置預(yù)設(shè)變形點(diǎn)對應(yīng)的幅值分布圖Fig.10 Amplitude distribution diagrams corresponding to preset deformation points at different positions

由圖10 仿真結(jié)果可知，當(dāng)點(diǎn)陣圖未發(fā)生點(diǎn)的偏移時(shí)，對圖像進(jìn)行傅里葉變換，濾波及傅里葉逆變換后，得到的像面幅值除點(diǎn)陣邊緣的突變產(chǎn)生的高頻信號外，點(diǎn)陣中央均勻排布的點(diǎn)陣經(jīng)一系列處理得到的像面幅值仍均勻排布，不存在峰值突變的地方，與前文分析的理論一致。而分別使點(diǎn)陣中間一點(diǎn)偏移、左上角一點(diǎn)偏移以及右下角一點(diǎn)偏移得到的點(diǎn)陣圖，經(jīng)相同的處理，得到的幅值圖中點(diǎn)偏移的部位會出現(xiàn)一個(gè)明顯的幅值突變，且幅值突變最大值處的坐標(biāo)剛好為原圖中點(diǎn)偏移后坐標(biāo)和未偏移的理論坐標(biāo)連線的中點(diǎn)，這說明最后得到的幅值分布圖可以直觀精確的定位發(fā)生點(diǎn)偏移的位置。

3.1.2 濾波頻點(diǎn)的選擇

分別選取左上角一點(diǎn)發(fā)生偏移的點(diǎn)陣圖得到的頻譜圖2 中的8 個(gè)頻點(diǎn)進(jìn)行濾波、逆傅里葉變換等操作得到其幅值分布圖，可觀察到選取關(guān)于原點(diǎn)對稱的一對頻點(diǎn)時(shí)得到的幅值分布圖完全相同。待測圖是左上角一點(diǎn)的偏移使得原本周期性分布的點(diǎn)陣在此處的空間頻率發(fā)生了改變，右上因?yàn)槠泣c(diǎn)的靠近而使得空間頻率更小，經(jīng)逆傅里葉變換得到的幅值圖中在該處應(yīng)有空間頻率發(fā)生變化的體現(xiàn)。此時(shí)選取高頻點(diǎn)1 和5 進(jìn)行濾波，得到的是一個(gè)包含波峰和波谷的幅值突變，右上方對應(yīng)了頻率更高的波峰，左下相應(yīng)的變?yōu)榈皖l的波谷；同樣的選取頻點(diǎn)2 和4 進(jìn)行濾波時(shí)得到的圖形相似，峰值有一個(gè)明顯向右向上的偏移（對應(yīng)于缺陷點(diǎn)位置的偏移），得到的峰值坐標(biāo)與理論值相比有1 到2 個(gè)像素的偏差；選取頻點(diǎn)3 進(jìn)行濾波時(shí)從圖中不能直觀反映缺陷點(diǎn)的偏移方向，但峰值基本居于正中心，更容易確定峰值點(diǎn)的坐標(biāo)，得到的坐標(biāo)與理論值幾乎沒有差別。因此，在已知只有一點(diǎn)發(fā)生偏移時(shí)，可通過選取合適的頻點(diǎn)來直觀地得到偏移的方向或者對偏移點(diǎn)進(jìn)行精確定位。

在探究單點(diǎn)偏移選取不同頻點(diǎn)對結(jié)果圖影響的基礎(chǔ)上，對預(yù)設(shè)同時(shí)包含有不同方向偏移點(diǎn)的點(diǎn)陣圖3 進(jìn)行傅里葉變換。圖11（a）、圖11（b）、圖11（c）、圖11（d）是分別選取頻點(diǎn)1、頻點(diǎn)2、頻點(diǎn)3、頻點(diǎn)4 進(jìn)行逆傅里葉變換得到的幅值分布圖。

圖11 不同方向偏移點(diǎn)選取不同頻點(diǎn)的幅值圖Fig.11 Amplitude diagrams of different frequency points selected at offset points in different directions

通過對比圖11（b）、圖11（c）、圖11（d）、圖11（e）幅圖可以發(fā)現(xiàn)，頻點(diǎn)2 只對水平方向的位移敏感，頻點(diǎn)4 只對豎直方向上的位移敏感，而對角線上的頻點(diǎn)1，3 對水平和豎直方向上的位移都敏感，不同位置的頻點(diǎn)對不同方向的位移敏感程度不同。因此，當(dāng)待測圖像中有一系列未知偏移方向的點(diǎn)時(shí)，為得到較為全面的表征一系列偏移點(diǎn)信息的幅值分布圖，應(yīng)對分別選取8 個(gè)高頻點(diǎn)得到的所有幅值圖相加取平均。

3.1.3 濾波器窗口大小

對左上角存在一點(diǎn)偏移的點(diǎn)陣圖進(jìn)行二維傅里葉變換，在得到的頻譜圖中選取同一個(gè)高頻點(diǎn)，只改變?yōu)V波時(shí)窗口的大小進(jìn)行濾波及逆傅里葉變換，記錄不同濾波窗口大小時(shí)從幅值分布圖中得到的偏移點(diǎn)的坐標(biāo)信息，再以原圖中偏移點(diǎn)偏移后坐標(biāo)和未偏移的理論坐標(biāo)連線中點(diǎn)處作為待測偏移點(diǎn)的理論坐標(biāo)值，繪制出不同濾波窗口大小下仿真實(shí)驗(yàn)得到的坐標(biāo)與理論坐標(biāo)的偏差，如圖12（b）所示。

圖12 濾波窗口大小對實(shí)驗(yàn)結(jié)果的影響Fig.12 Influence of filter window size on experimental results

顯然，在濾波窗口取60 像素到85 像素時(shí)，幾乎不存在誤差；濾波窗口取45 像素到55 像素時(shí)也僅存在1 像素到2 像素的誤差；濾波窗口更小時(shí)（對應(yīng)于頻譜圖中的A 圈），由于高頻點(diǎn)的信息丟失的較多因而誤差較大；同時(shí)，當(dāng)濾波窗口取到90 像素時(shí)（對應(yīng)于頻譜圖中的B 圈），由于窗口太大引入了相鄰高頻點(diǎn)的信息，對原信息造成了干擾因而出現(xiàn)了較大誤差。

仿真發(fā)現(xiàn)，改變?yōu)V波器窗口的大小，得到的幅值分布圖在細(xì)節(jié)體現(xiàn)上會存在差異，但只要保證濾波器的窗口能包含高頻點(diǎn)的全部信息且沒有將相鄰高頻點(diǎn)的信息混入，對同一個(gè)高頻點(diǎn)進(jìn)行濾波，再做逆傅里葉變換，得到的幅值分布圖中表征偏移點(diǎn)位置的幅值突變處的坐標(biāo)差別并不大。

3.1.4 偏移點(diǎn)的大小

對不同大小但圓心坐標(biāo)相同的偏移點(diǎn)經(jīng)過傅里葉變換，選取同一高頻點(diǎn)以同一窗口大小對其進(jìn)行濾波，再逆傅里葉變換得到幅值分布圖，從得到的幅值分布圖中提取偏移點(diǎn)的位置信息并與理論值作比較分別繪出圖13（a）、圖13（b）、圖13（c）。

圖13 偏移點(diǎn)大小對實(shí)驗(yàn)結(jié)果的影響Fig.13 Influence of offset point size on experimental results

圖13（a）和圖13（b）圖分別為偏移點(diǎn)取不同直徑時(shí)對應(yīng)的x軸、y軸的偏差，可以看出，當(dāng)偏移點(diǎn)的直徑取15 像素時(shí)（即偏移點(diǎn)與其余未偏移點(diǎn)大小一樣），得到的偏移點(diǎn)的位置信息與理論值完全相同，在偏移點(diǎn)直徑進(jìn)一步減小到13 像素時(shí)，x軸、y軸方向仍沒有出現(xiàn)偏差，偏移點(diǎn)直徑減小到11 像素時(shí)也僅有1 像素的偏差。這說明當(dāng)偏移點(diǎn)直徑減少的值在標(biāo)準(zhǔn)偏移點(diǎn)直徑的26.7%的范圍內(nèi)，對測量結(jié)果不會產(chǎn)生太大的影響，該測量方法對于偏移點(diǎn)減小的情況有一定的容錯(cuò)率。但是當(dāng)偏移點(diǎn)進(jìn)一步減小時(shí)，x軸、y軸都開始出現(xiàn)偏差，且偏移點(diǎn)越小兩軸出現(xiàn)的偏差越大。通過圖13（a）和圖13（b）兩圖還可以發(fā)現(xiàn)，隨著偏移點(diǎn)直徑減小，x軸方向的偏差減小到負(fù)值，而y軸方向的偏差增大到正值，結(jié)合圖偏移點(diǎn)向右向上的偏移方向可以看出，當(dāng)偏移點(diǎn)進(jìn)一步減小時(shí)，測量得到的坐標(biāo)值與理論坐標(biāo)值相比，產(chǎn)生了一個(gè)向左、向下的位置偏差，正好與偏移點(diǎn)偏移方向相反。而偏移點(diǎn)直徑增大時(shí)所有的x軸方向的偏差、y軸方向的偏差較小且均為定值。綜上，當(dāng)偏移點(diǎn)的大小在標(biāo)準(zhǔn)直徑的26.7%的范圍內(nèi)減小時(shí)，對測量結(jié)果不會產(chǎn)生太大的影響，減小值超出這一范圍時(shí)，測得坐標(biāo)值會削弱偏移點(diǎn)產(chǎn)生的位置偏移量，使測得的位置偏移量比實(shí)際的小。移點(diǎn)當(dāng)偏移點(diǎn)增大此時(shí)的偏差幾乎可以忽略且不隨偏移點(diǎn)直徑的改變而改變。實(shí)際中微柱測細(xì)胞牽引力時(shí)對應(yīng)于偏移點(diǎn)減小的情況，只要滿足微柱頂端投影的減少量在微柱標(biāo)準(zhǔn)直徑的26.7%范圍內(nèi)，即使被拉伸的微柱頂端投影相較于其他未偏移的微柱頂端直徑小，仍能得到比較準(zhǔn)確的位置信息；微柱頂端投影直徑的減少量超過微柱標(biāo)準(zhǔn)直徑的26.7%時(shí)會使測得的微柱頂端偏移量偏小。

3.2 細(xì)胞-微柱顯微圖

圖14 為經(jīng)填充的微柱高倍顯微灰度圖進(jìn)行一次傅里葉變換得到的頻譜圖。分別選取8 個(gè)一級衍射斑點(diǎn)進(jìn)行傅里葉逆變換，將得到的8 張幅值分布圖相加取平均作為最終的幅值圖。通過偏移值與幅值的對應(yīng)關(guān)系將幅值分布轉(zhuǎn)換成微柱偏移量的分布圖，再由PDMS 的剛度計(jì)算得到細(xì)胞力的分布圖如圖15（a）所示。

圖14 進(jìn)行傅里葉變換得到的頻譜圖Fig.14 Obtained spectrogram by Fourier transform

圖15 細(xì)胞顯微圖的處理Fig.15 Cell micrograph processing

對同一幅高倍顯微圖用質(zhì)心法測細(xì)胞力，先將其轉(zhuǎn)換為灰度圖，再進(jìn)行中值濾波、二值化等處理，利用Matlab 中regionprops 函數(shù)提微柱端面質(zhì)心坐標(biāo)。將沒有細(xì)胞粘附處未發(fā)生偏轉(zhuǎn)的微柱坐標(biāo)作為參考坐標(biāo)，粘附有細(xì)胞的偏轉(zhuǎn)微柱頂端質(zhì)心坐標(biāo)與參考坐標(biāo)作差得到微柱的偏移量分布圖，再乘仿真得到的剛度得到細(xì)胞力分布如圖15（b）。

對比圖15 兩幅圖可以看出，傅里葉變換的方法與質(zhì)心法得到細(xì)胞力的分布圖基本一致，從兩種方法得到的細(xì)胞力分布圖中各選取被細(xì)胞覆蓋且通過牽引力帶動偏轉(zhuǎn)的部分，并將其分為5 個(gè)相對應(yīng)的主要區(qū)域進(jìn)行研究，在圖15 中對同一區(qū)域進(jìn)行標(biāo)號。以質(zhì)心法得到的結(jié)果作為標(biāo)準(zhǔn)值，分別計(jì)算出傅里葉變換法得到的5 個(gè)區(qū)域的最大細(xì)胞力，如表1 所示。

表1 測得的各區(qū)域最大細(xì)胞力Table 1 Maximum cell forces measured in each region

從表1 可以看出，傅里葉變換法得到的細(xì)胞力分布圖中各區(qū)域細(xì)胞力的最大值與質(zhì)心法相比，相對誤差在1.76%到17.37%的范圍內(nèi)，且總體測得的細(xì)胞力均小于質(zhì)心法測得的細(xì)胞力。這種情況的產(chǎn)生可能是由于傅里葉變換法幅值分布圖是由分別取各個(gè)高頻點(diǎn)得到的幅值圖相加取平均得到的，而如前文分析，某些頻點(diǎn)只對特定方向的偏移敏感，這樣得到的頻點(diǎn)的幅值分布圖明顯會削弱其他方向的偏移信息，因此所有頻點(diǎn)相加取平均，測得的總體值會偏小。此外，區(qū)域3、區(qū)域4 的相對誤差分別為16.54%和17.37%，相對其他區(qū)域的相對誤差值較大，測得的細(xì)胞力較質(zhì)心法測得的細(xì)胞力偏小的程度更明顯。這種情況可能是這兩個(gè)區(qū)域細(xì)胞牽引力較大而使微柱發(fā)生了相對較大的位移，從而使得微柱頂端的投影直徑相比較于其他微柱直徑減小量較多造成的。從圖5 細(xì)胞的高倍顯微圖中也可以看出，此時(shí)這兩個(gè)區(qū)域的微柱頂端均發(fā)生了較大的位移，微柱頂端面積約為其他微柱頂端面積的一半。而如前文分析，偏移點(diǎn)直徑的減小量超過標(biāo)準(zhǔn)偏移點(diǎn)直徑的26.7%，導(dǎo)致了測量結(jié)果偏小。

表2 中，傅里葉變換法得到的細(xì)胞力分布圖對各區(qū)域的細(xì)胞力取平均值，各區(qū)域相對誤差值與表1 相比方差更小，測量的精度受偏移點(diǎn)偏移方向的影響更小，該方法測得的各區(qū)域的平均細(xì)胞力誤差更小。

表2 測得的各區(qū)平均細(xì)胞力Table 2 Mean cell forces measured in each region

此外，分別用數(shù)字傅里葉變換和質(zhì)心法對處理點(diǎn)陣個(gè)數(shù)分別為900（30×30）、3600（60×60）、8100（90×90）二維點(diǎn)陣圖的運(yùn)行時(shí)間如圖15所示。

從圖中16 可以明顯看到，在二維點(diǎn)陣圖中點(diǎn)的個(gè)數(shù)較小時(shí)，質(zhì)心法的運(yùn)行時(shí)間和傅里葉變換相差不大；當(dāng)待處理的點(diǎn)數(shù)繼續(xù)增大（從30×30 的點(diǎn)陣增大到60×60 的點(diǎn)陣）時(shí)，傅里葉變換的運(yùn)行時(shí)間基本不變，但質(zhì)心法的運(yùn)行時(shí)間增加近2 倍；待處理點(diǎn)數(shù)增加到8100（90×90）時(shí)，傅里葉變換運(yùn)行時(shí)間仍然基本不變，但是質(zhì)心法的運(yùn)行時(shí)間增加到了數(shù)字傅里葉變換的10 倍。由此可以發(fā)現(xiàn)，在處理待測點(diǎn)數(shù)較多的細(xì)胞顯微圖時(shí)，傅里葉變換法相較于質(zhì)心法，運(yùn)行時(shí)間更短，從執(zhí)行速度上來說該方法也具有明顯優(yōu)勢。

圖16 傅里葉變和質(zhì)心法運(yùn)行時(shí)間的比較Fig.16 Comparison of running time between Fourier transform and centroid method

4 結(jié)論

本文提出了一種基于傅里葉變換的細(xì)胞牽引力測量方法，并通過對預(yù)設(shè)二維點(diǎn)陣的仿真驗(yàn)證了測量方法的可行性，并對過程中頻點(diǎn)選擇、濾波窗口大小進(jìn)行分析，確定了最優(yōu)值。在此基礎(chǔ)上，分別用傅里葉變換法和質(zhì)心法對同一細(xì)胞-微柱顯微圖進(jìn)行了測量，傅里葉變換法的細(xì)胞力分布圖在結(jié)果呈現(xiàn)上不亞于質(zhì)心法，測得的各區(qū)域最大細(xì)胞力相對誤差在17.37%之內(nèi)，各區(qū)域平均細(xì)胞力相對誤差在7.93%之內(nèi)。該方法避免了照明不均時(shí)閾值的合理選取困難、需原始位置作為標(biāo)準(zhǔn)、二值化的過程使得輪廓跳動大這些質(zhì)心法中存在的問題，且對于8100 個(gè)點(diǎn)數(shù)的點(diǎn)陣圖的運(yùn)算速度比質(zhì)心法快近10 倍，在運(yùn)算速度上具有明顯優(yōu)勢。此外，該方法可以直觀體現(xiàn)細(xì)胞力的總體分布情況，為后續(xù)通過光學(xué)傅里葉變換實(shí)現(xiàn)細(xì)胞力的測量提供了研究基礎(chǔ)，對光學(xué)傅里葉變換在大視場下同時(shí)對多個(gè)細(xì)胞的細(xì)胞牽引力分布情況進(jìn)行研究有重要的參考意義。