專利名稱：一種細胞原位主動變形測量方法
技術(shù)領域：
本發(fā)明涉及一種細胞顯微變形測量方法，特別是關(guān)于一種用于定量分析細胞與細胞外基質(zhì)(彈性基底)之間物理相互作用的細胞原位主動變形測量方法。

背景技術(shù)：
細胞是研究生命科學的基礎，也是現(xiàn)代生命科學發(fā)展的重要支柱。一般來說，動物組織中的細胞不僅與相鄰細胞接觸和作用，同時也與由分泌蛋白和多糖組成的充滿胞外空間的復雜網(wǎng)絡結(jié)構(gòu)，即細胞外基質(zhì)相互接觸和作用。細胞與細胞外基質(zhì)等物理環(huán)境相互作用涉及許多復雜的生化信號轉(zhuǎn)導通路，定量研究這種相互作用行為具有重要的生理病理意義。研究細胞與胞外基質(zhì)相互作用的一類常用方法是采用彈性硅膠基底膜來模擬細胞外基質(zhì)微環(huán)境，把細胞培養(yǎng)在其上面，通過牽拉硅膠膜的方式促使細胞發(fā)生被動變形，以此研究胞外基質(zhì)環(huán)境變化對與細胞粘附、鋪展及生長、分化等生物學行為的影響。1999年，Dembo(丹博)和王毓立等人通過把細胞培養(yǎng)在表面摻有熒光顆粒的聚丙烯酰胺薄膜彈性基底上，觀察到了由細胞主動變形(包括細胞的周期性收縮和遷移)所導致的基底膜微納尺度形變，以此為基礎研究細胞的自主收縮、定向遷移規(guī)律，探索單細胞對外界環(huán)境的感知機理(Dembo M，Wang YL.Stresses at the cell-to-substrate interface duringlocomotion of fibroblasts.Biophysical Journal，1999；762307-2316.丹博，王毓立，成纖維細胞遷移過程中細胞-基底界面應力表征，生物物理，1999，762307-2316)。
目前，定量測量這種薄膜變形的方法主要有基于光流的區(qū)域特征匹配方法(Marganski WA，Dembo M，Wang YL.Measurements of cell-generated deformationson flexible substrata using correlation-based optical flow.Methods Enzymol.2003；361197-211.馬甘斯基，丹博，王毓立，運用基于相關(guān)的光流方法測量細胞導致的軟基底變形，酶學方法，2003；361197-211.)以及基于模式識別的熒光粒子點跟蹤方法(Yang ZC，Lin JS，Chen JX，Wang J.Determining substrate displacement andcell traction fields-a new approach.J.Theor.Biol.2006；242607-616.楊兆春，林津尚，陳建新，王慧聰，確定細胞牽引力場——一種新的方法，理論生物學，2006；242607-616)。然而，基于光流的區(qū)域特征匹配方法搜索計算量大，計算速度很難滿足要求；基于模式識別的熒光粒子點跟蹤方法計算速度稍快，但由于單個的熒光粒子光斑幾何光學特征相同，實際處理過程中極易出現(xiàn)錯誤跟蹤，計算結(jié)果的準確性很難保證。


發(fā)明內(nèi)容
針對上述問題，本發(fā)明的目的是提供一種具有較高的準確率、效率及高時間-空間分辨率的細胞原位主動變形測量方法。
為實現(xiàn)上述目的，本發(fā)明采取以下技術(shù)方案一種細胞原位主動變形測量方法，其包括以下步驟(1)分別在未發(fā)生彈性變形及發(fā)生彈性變形的細胞基底熒光圖的參考圖像F(x，y)和變形圖像G(x，y)上，以設定的采樣間距S為邊長劃分方形采樣網(wǎng)格，單位為像素；(2)分別以各采樣點(α，β)為中心，在參考圖像F(x，y)上設定邊長為(2N+1)的方形相關(guān)子區(qū)fα，β(x，y)，在變形圖像G(x，y)上選擇邊長為(2M+1)的方形搜索區(qū)域Gα，β(x，y)，其中N的取值范圍為15～30之間的整數(shù)，M的取值范圍為20～80之間的整數(shù)，單位均為像素，且M＞N；(3)根據(jù)參考圖像F(x，y)和變形圖像G(x，y)的數(shù)字圖像灰度分布結(jié)構(gòu)表f(x，y)和g(x，y)，利用遞推方法分別建立圖像灰度全局求和結(jié)構(gòu)表Sf、Sg及圖像能量全局求和結(jié)構(gòu)表


(4)根據(jù)數(shù)字圖像灰度分布結(jié)構(gòu)表f(x，y)和g(x，y)，計算以采樣點(α，β)為中心的相關(guān)子區(qū)fα，β(x，y)與相應的搜索區(qū)域Gα，β(x，y)之間的互相關(guān)系數(shù)矩陣Pα，β，得到其快速遞推關(guān)系式為Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)，其中，令

(u，v)表示參考圖像的相關(guān)子區(qū)fα，β(x，y)中采樣點(α，β)在變形圖像的搜索區(qū)域Gα，β(x，y)中的位移；D是相關(guān)子區(qū)fα，β(x，y)的邊長，D＝2N+1；(5)根據(jù)所述步驟(3)、(4)的結(jié)果，計算參考圖像的相關(guān)子區(qū)fα，β(x，y)與相應變形圖像的搜索區(qū)域Gα，β(x，y)的變形子區(qū)gα，β(x，y)之間的零均值歸一化互相關(guān)系數(shù)矩陣Gα，β(u，v)；(6)采用基于梯度的亞像素位移定位算法，在零均值歸一化互相關(guān)系數(shù)矩陣的最大峰值所在位置附近進行亞像素插值運算，求得參考圖像中采樣點(α，β)在變形圖像中的準確位置為

式中，(Uα，β，Vα，β)表示采樣點(α，β)在變形圖像中的總位移大小，(uα，β，vα，β)是由所述零均值歸一化相關(guān)系數(shù)矩陣Cα，β(u，v)中最大元素所確定的整像素位移，(Δuα，β，Δvα，β)是(uα，β，vα，β)相應的亞像素位移；(7)重復所述步驟(4)～步驟(6)，計算所有采樣點在變形圖像中的準確位置，進而得到整張變形圖像位移場。
所述步驟(4)中，所述互相關(guān)系數(shù)矩陣Pα，β快速遞推關(guān)系式由以下步驟得到①對于采樣點(α，β)，令相關(guān)子區(qū)fα，β(x，y)與相應的搜索區(qū)域Gα，β(x，y)之間構(gòu)成互相關(guān)系數(shù)矩陣為Pα，β；②與采樣點(α，β)相鄰的采樣點(α+S，β)，其對應的相關(guān)子區(qū)fα+S，β(x，y)與相應的搜索區(qū)域Gα+S，β(x，y)之間的互相關(guān)系數(shù)矩陣為Pα+S，β；③將兩個互相關(guān)系數(shù)矩陣Pα，β和Pα+S，β對比得到Pα+S，β(u，v)＝Pα，β(u，v)-Jα，β(u，v)+Iα，β(u，v)，其中，令

得到Jα，β(u，v)＝Iα-D，β(u，v)，D是相關(guān)子區(qū)fα，β(x，y)的邊長，D＝2N+1；④由所述步驟①～③中各式能建立快速遞推關(guān)系為Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)。
所述步驟(5)中，所述零均值歸一化互相關(guān)系數(shù)矩陣Cα，β(u，v)為令其中，P(α，β；u，v)為圖像參考子區(qū)fα，β(x，y)與相應搜索區(qū)域Gα，β(x，y)之間的互相關(guān)系數(shù)矩陣，其由所述步驟(4)中的所述互相關(guān)系數(shù)矩陣Pα，β快速遞推關(guān)系式得到；F(α，β)為參考圖像灰度平方和；G(α，β；u，v)為變形圖像灰度平方和；Q(α，β；u，v)、F(α，β)和G(α，β；u，v)通過快速查詢所述步驟(3)中預先建立的四個圖像灰度、圖像能量全局求和結(jié)構(gòu)表得到。
本發(fā)明由于采取以上技術(shù)方案，其具有以下優(yōu)點1、本發(fā)明由于采用了基于快速遞推關(guān)系式和快速查詢四個圖像灰度、圖像能量全局求和結(jié)構(gòu)表的方式，快速計算零均值歸一化相關(guān)系數(shù)矩陣，并根據(jù)該矩陣的峰值位置確定采樣點的整像素位移大小，進而得到采樣點的位移，因此實現(xiàn)了測量精確高效、高時間-空間分辨率的效果。2、本發(fā)明由于采用了快速查表方式和快速遞推方法進行計算，因此可以有效地降低采樣點位移計算的復雜性，特別是在高時間空間位移采樣以及較大搜索區(qū)域尺寸情況下，能顯著的提高細胞-基底位移場求解效率。3、本發(fā)明由于采用了快速遞推方式進行計算，這種方法比傳統(tǒng)數(shù)字圖像相關(guān)方法的計算效率大約提高了10到50倍，特別是當采樣網(wǎng)格邊長減小時，計算效率尤為突出。因此，本發(fā)明可以應用于高時空分辨率的細胞原位主動變形測量中。



圖1是本發(fā)明的整體流程示意圖 圖2是本發(fā)明采樣點網(wǎng)格劃分的相關(guān)子區(qū)、搜索區(qū)域示意圖 圖3是本發(fā)明的實施例一中單細胞粘附在聚丙烯酰胺彈性基底上的相差圖 圖4是本發(fā)明的實施例一中聚丙烯酰胺基底的變形前參考圖像 圖5是本發(fā)明的實施例一中聚丙烯酰胺基底的變形圖像 圖6是本發(fā)明的實施例一中基于本發(fā)明測量方法計算得到的熒光基底位移場
具體實施例方式 下面結(jié)合附圖和實施例對本發(fā)明進行詳細的描述。
如圖1所示，本發(fā)明是基于零均值歸一化互相關(guān)系數(shù)在細胞基底位移場求解過程中所表現(xiàn)出的內(nèi)在特點，通過構(gòu)建與圖像灰度和圖像能量密切相關(guān)的四個全局求和表，發(fā)展了一整套位移場求解快速遞推方法，其步驟如下 1)將未發(fā)生彈性變形的細胞基底熒光圖作為參考圖像F(x，y)，發(fā)生彈性變形的細胞基底熒光圖作為變形圖像G(x，y)；設定采樣間距為S，分別在參考圖像F(x，y)和變形圖像G(x，y)上以S為邊長劃分規(guī)則的方形采樣網(wǎng)格(如圖2所示)，單位為像素； 2)分別以每個采樣點(α，β)為中心，在參考圖像F(x，y)上設定邊長為(2N+1)的方形相關(guān)子區(qū)fα，β(x，y)，在變形圖像G(x，y)上選擇邊長為(2M+1)的方形搜索區(qū)域Gα，β(x，y)，其中N的取值范圍一般為15～30之間的任意整數(shù)(其單位為像素)，而M的取值范圍一般為20～80之間的任意整數(shù)(其單位也為像素)，且M＞N； 3)根據(jù)參考圖像F(x，y)和變形圖像G(x，y)的數(shù)字圖像灰度分布結(jié)構(gòu)表f(x，y)和g(x，y)(如表1、表2所示)，利用遞推方法分別建立參考圖像F(x，y)和變形圖像G(x，y)的圖像灰度全局求和結(jié)構(gòu)表Sf、Sg及圖像能量全局求和結(jié)構(gòu)表


如表3～表6所示，其遞推方式如下 Sf(x，y)＝f(x，y)+Sf(x-1，y)+Sf(x，y-1)-Sf(x-1，y-1)，(1) Sg(x，y)＝g(x，y)+Sg(x-1，y)+Sg(x，y-1)-Sg(x-1，y-1)， (3) 其中，當x，y≤0時，
表1 參考數(shù)字圖像灰度分布結(jié)構(gòu)表f(x，y) 表2 變形后數(shù)字圖像灰度分布結(jié)構(gòu)表g(x，y) 表3 參考圖相灰度求和結(jié)構(gòu)表Sf 表4 變形圖相灰度求和結(jié)構(gòu)表Sg 表5 參考圖像能量求和結(jié)構(gòu)表


表6 參考圖像能量求和結(jié)構(gòu)表

上述各表格中，p、q、w、H均為參考的圖像參考子區(qū)fα，β(x，y)與相應搜索區(qū)域Gα，β(x，y)內(nèi)任意點的取值，其單位為像素； 4)運用遞推方法，根據(jù)參考圖像F(x，y)和變形圖像G(x，y)的數(shù)字圖像灰度分布結(jié)構(gòu)表f(x，y)和g(x，y)，計算以采樣點(α，β)為中心的參考圖像的相關(guān)子區(qū)fα，β(x，y)與相應變形圖像的搜索區(qū)域Gα，β(x，y)之間的互相關(guān)系數(shù)矩陣Pα，β，進而得到互相關(guān)系數(shù)矩陣Pα，β快速遞推關(guān)系式為 Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)，(5) 令 其中，(u，v)表示參考圖像的相關(guān)子區(qū)fα，β(x，y)中采樣點(α，β)在變形圖像的搜索區(qū)域Gα，β(x，y)中的位移(如圖2所示)；D是相關(guān)子區(qū)fα，β(x，y)的邊長，D＝2N+1； 5)根據(jù)步驟3)和步驟4)的結(jié)果，計算參考圖像的相關(guān)子區(qū)fα，β(x，y)與相應變形圖像的搜索區(qū)域Gα，β(x，y)的變形子區(qū)gα，β(x，y)之間的零均值歸一化互相關(guān)系數(shù)矩陣(Zero-normalized cross-correlation coefficient，ZNCC)，該零均值歸一化互相關(guān)系數(shù)矩陣Cα，β(u，v)為 令 其中，P(α，β；u，v)為圖像參考子區(qū)fα，β(x，y)與相應搜索區(qū)域Gα，β(x，y)之間的互相關(guān)系數(shù)矩陣；F(α，β)為參考圖像灰度平方和；G(α，β；u，v)為變形圖像灰度平方和； 公式(7)可由標準的零均值歸一化互相關(guān)系數(shù)表達式經(jīng)化簡得到下式 上述公式中，

和

分別表示參考圖像灰度和變形圖像灰度的均值； 6)采用基于梯度的亞像素位移定位算法，在步驟5)中的ZNCC系數(shù)矩陣最大峰值所在位置附近進行亞像素插值運算，這樣可以把熒光彈性基底的位移測量精度由整像素水平延伸至亞像素水平，即求得參考圖像中采樣點(α，β)在變形圖像中的準確位置為 上式中，(Uα，β，Vα，β)表示采樣點(α，β)在變形圖像中的總位移大小，(uα，β，vα，β)是由步驟5)中的零均值歸一化相關(guān)系數(shù)矩陣Cα，β(u，v)中最大元素所確定的整像素位移，(Δuα，β，Δvα，β)是(uα，β，vα，β)相應的亞像素位移； 7)重復步驟4)～步驟6)，計算所有采樣點在變形后圖像中的準確位置，進而得到整張變形圖像位移場。
上述步驟4)中，互相關(guān)系數(shù)矩陣Pα，β快速遞推關(guān)系式通過以下步驟得到 ①對于采樣點(α，β)，令參考圖像的相關(guān)子區(qū)fα，β(x，y)與相應的變形圖像的搜索區(qū)域Gα，β(x，y)之間構(gòu)成的互相關(guān)系數(shù)矩陣Pα，β為
②同理，對于與采樣點(α，β)相鄰的采樣點(α+S，β)而言，其對應的相關(guān)子區(qū)fα+S，β(x，y)與相應的搜索區(qū)域Gα+S，β(x，y)之間的互相關(guān)系數(shù)矩陣Pα+S，β可以表示為
③將上述兩個互相關(guān)系數(shù)矩陣Pα，β和Pα+S，β進行對比后，可以得到 其中，令 顯然，將Jα，β(u，v)與公式(6)進行對比后可知，Jα，β(u，v)＝Iα-D，β(u，v)，D是相關(guān)子區(qū)fα，β(x，y)的邊長，即D＝2N+1； ④由上述步驟中各式可建立如下快速遞推關(guān)系 Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)。
上述步驟5)中，互相關(guān)系數(shù)矩陣P(α，β；u，v)根據(jù)步驟4)中的互相關(guān)系數(shù)矩陣Pα，β快速遞推關(guān)系式可以得到 公式(7)中的Q(α，β；u，v)、F(α，β)和G(α，β；u，v)可以通過快速查詢步驟3)中預先建立的四個圖像灰度、圖像能量全局求和結(jié)構(gòu)表(表3～表6)得到，即 上述步驟6)中，亞像素位移(Δuα，β，Δvα，β)為 上式中，

和

分別表示在u和v方向上灰度的一階梯度，灰度梯度算子取為[1/12，-8/12，0，8/12，-1/12]，這時 下面通過具體實施例對本發(fā)明的測量方法進行進一步的描述。
實施例一如圖3～圖6所示，給出的一組細胞變形圖像，其圖像尺寸為1027×533像素，當參考子區(qū)邊長為40和60個像素時的計算效率測試結(jié)果如表7、表8所示， 表7 傳統(tǒng)數(shù)字圖像相關(guān)方法與本發(fā)明的測量方法計算時間之比(即加速比)
表8 傳統(tǒng)數(shù)字圖像相關(guān)方法與本發(fā)明的測量方法計算時間之比(即加速比)
從上述兩個表中可以看到，當搜索半徑介于10～50個像素，采樣網(wǎng)格邊長為5～20個像素時，傳統(tǒng)數(shù)字圖像相關(guān)方法與本發(fā)明的測量方法計算時間之比(即加速比)大約在10～50之間變化，這表明本發(fā)明的測量方法能夠較傳統(tǒng)數(shù)字圖像相關(guān)方法提高計算效率大約10到50倍。特別是當采樣網(wǎng)格邊長減小時，計算效率尤為突出，這一特點預示著本發(fā)明的測量方法更適于高時空分辨率的細胞原位主動變形測量。
上述各實施例僅用于說明本發(fā)明，凡是在本發(fā)明技術(shù)方案的基礎上進行的等同變換和改進，均不應排除在本發(fā)明的保護范圍之外。
權(quán)利要求
1.一種細胞原位主動變形測量方法，其包括以下步驟
(1)分別在未發(fā)生彈性變形及發(fā)生彈性變形的細胞基底熒光圖的參考圖像F(x，y)和變形圖像G(x，y)上，以設定的采樣間距S為邊長劃分方形采樣網(wǎng)格，單位為像素；
(2)分別以各采樣點(α，β)為中心，在參考圖像F(x，y)上設定邊長為(2N+1)的方形相關(guān)子區(qū)fα，β(x，y)，在變形圖像G(x，y)上選擇邊長為(2M+1)的方形搜索區(qū)域Gα，β(x，y)，其中N的取值范圍為15～30之間的整數(shù)，M的取值范圍為20～80之間的整數(shù)，單位均為像素，且M＞N；
(3)根據(jù)參考圖像F(x，y)和變形圖像G(x，y)的數(shù)字圖像灰度分布結(jié)構(gòu)表f(x，y)和g(x，y)，利用遞推方法分別建立圖像灰度全局求和結(jié)構(gòu)表Sf、Sg及圖像能量全局求和結(jié)構(gòu)表
(4)根據(jù)數(shù)字圖像灰度分布結(jié)構(gòu)表f(x，y)和g(x，y)，計算以采樣點(α，β)為中心的參考圖像的相關(guān)子區(qū)fα，β(x，y)與相應變形圖像的搜索區(qū)域Gα，β(x，y)之間的互相關(guān)系數(shù)矩陣Pα，β，得到其快速遞推關(guān)系式為
Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)，
其中，令
(u，v)表示參考圖像的相關(guān)子區(qū)fα，β(x，y)中采樣點(α，β)在變形圖像的搜索區(qū)域Gα，β(x，y)中的位移；D是相關(guān)子區(qū)fα，β(x，y)的邊長，D＝2N+1；
(5)根據(jù)所述步驟(3)、(4)的結(jié)果，計算參考圖像的相關(guān)子區(qū)fα，β(x，y)與相應變形圖像的搜索區(qū)域Gα，β(x，y)的變形子區(qū)gα，β(x，y)之間的零均值歸一化互相關(guān)系數(shù)矩陣Cα，β(u，v)；
(6)采用基于梯度的亞像素位移定位算法，在零均值歸一化互相關(guān)系數(shù)矩陣的最大峰值所在位置附近進行亞像素插值運算，求得參考圖像中采樣點(α，β)在變形圖像中的準確位置為
式中，(Uα，β，Vα，β)表示采樣點(α，β)在變形圖像中的總位移大小，(uα，β，vα，β)是由所述零均值歸一化相關(guān)系數(shù)矩陣Cα，β(u，v)中最大元素所確定的整像素位移，(Δuα，β，Δvα，β)是(uα，β，vα，β)相應的亞像素位移；
(7)重復所述步驟(4)～步驟(6)，計算所有采樣點在變形圖像中的準確位置，得到整張變形圖像位移場。
2.如權(quán)利要求1所述的一種細胞原位主動變形測量方法，其特征在于所述步驟(4)中，所述互相關(guān)系數(shù)矩陣Pα，β快速遞推關(guān)系式由以下步驟得到
①對于采樣點(α，β)，令相關(guān)子區(qū)fα，β(x，y)與相應的搜索區(qū)域Gα，β(x，y)之間構(gòu)成互相關(guān)系數(shù)矩陣為Pα，β；
②與采樣點(α，β)相鄰的采樣點(α+S，β)，其對應的相關(guān)子區(qū)fα+S，β(x，y)與相應的搜索區(qū)域Gα+S，β(x，y)之間的互相關(guān)系數(shù)矩陣為Pα+S，β；
③將所述兩個互相關(guān)系數(shù)矩陣Pα，β和Pα+S，β進行對比后，得到
Pα+S，β(u，v)＝Pα，β(u，v)-Jα，β(u，v)+Iα，β(u，v)
其中，令
D是相關(guān)子區(qū)fα，β(x，y)的邊長，D＝2N+1；
④由所述步驟①～③中各式能建立快速遞推關(guān)系為
Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)。
3.如權(quán)利要求1所述的一種細胞原位主動變形測量方法，其特征在于所述步驟(5)中，所述零均值歸一化互相關(guān)系數(shù)矩陣Cα，β(u，v)為
令
其中，P(α，β；u，v)為圖像參考子區(qū)fα，β(x，y)與相應搜索區(qū)域Gα，β(x，y)之間的互相關(guān)系數(shù)矩陣，其由所述步驟(4)中的所述互相關(guān)系數(shù)矩陣Pα，β快速遞推關(guān)系式得到；F(α，β)為參考圖像灰度平方和；G(α，β；u，v)為變形圖像灰度平方和；Q(α，β；u，v)、F(α，β)和G(α，β；u，v)通過快速查詢所述步驟(3)中預先建立的四個圖像灰度、圖像能量全局求和結(jié)構(gòu)表得到。
全文摘要
本發(fā)明涉及一種細胞原位主動變形測量方法，其測量步驟包括(1)分別在細胞-基底熒光參考圖像和相應的變形圖像上劃分方形采樣網(wǎng)格；(2)確定參考子區(qū)和搜索區(qū)域的尺寸；(3)建立四個全局求和表；(4)以遞推方式計算以采樣點為中心的圖像參考子區(qū)與相應搜索區(qū)域之間的互相關(guān)系數(shù)；(5)利用步驟(4)的結(jié)果并通過查求和表的方式快速計算零均值歸一化相關(guān)系數(shù)矩陣，并根據(jù)此矩陣的峰值位置確定采樣點的整像素位移大??；(6)利用梯度法進行亞像素插值運算；(7)遞推計算所有采樣點處的位移，得到整張變形圖像精確位移場。本發(fā)明具有測量精確高效、高時間-空間分辨率的效果。本發(fā)明可以應用于高時空分辨率的細胞原位主動變形測量中。
文檔編號G01N33/48GK101813693SQ20101017069
公開日2010年8月25日 申請日期2010年5月6日 優(yōu)先權(quán)日2010年5月6日
發(fā)明者黃建永, 潘曉暢, 秦雷, 朱濤, 熊春陽, 方競 申請人:北京大學