专利名称：一种细胞原位主动变形测量方法
技术领域：
本发明涉及一种细胞显微变形测量方法，特别是关于一种用于定量分析细胞与细胞外基质(弹性基底)之间物理相互作用的细胞原位主动变形测量方法。

背景技术：
细胞是研究生命科学的基础，也是现代生命科学发展的重要支柱。一般来说，动物组织中的细胞不仅与相邻细胞接触和作用，同时也与由分泌蛋白和多糖组成的充满胞外空间的复杂网络结构，即细胞外基质相互接触和作用。细胞与细胞外基质等物理环境相互作用涉及许多复杂的生化信号转导通路，定量研究这种相互作用行为具有重要的生理病理意义。研究细胞与胞外基质相互作用的一类常用方法是采用弹性硅胶基底膜来模拟细胞外基质微环境，把细胞培养在其上面，通过牵拉硅胶膜的方式促使细胞发生被动变形，以此研究胞外基质环境变化对与细胞粘附、铺展及生长、分化等生物学行为的影响。1999年，Dembo(丹博)和王毓立等人通过把细胞培养在表面掺有荧光颗粒的聚丙烯酰胺薄膜弹性基底上，观察到了由细胞主动变形(包括细胞的周期性收缩和迁移)所导致的基底膜微纳尺度形变，以此为基础研究细胞的自主收缩、定向迁移规律，探索单细胞对外界环境的感知机理(Dembo M，Wang YL.Stresses at the cell-to-substrate interface duringlocomotion of fibroblasts.Biophysical Journal，1999；762307-2316.丹博，王毓立，成纤维细胞迁移过程中细胞-基底界面应力表征，生物物理，1999，762307-2316)。
目前，定量测量这种薄膜变形的方法主要有基于光流的区域特征匹配方法(Marganski WA，Dembo M，Wang YL.Measurements of cell-generated deformationson flexible substrata using correlation-based optical flow.Methods Enzymol.2003；361197-211.马甘斯基，丹博，王毓立，运用基于相关的光流方法测量细胞导致的软基底变形，酶学方法，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)。然而，基于光流的区域特征匹配方法搜索计算量大，计算速度很难满足要求；基于模式识别的荧光粒子点跟踪方法计算速度稍快，但由于单个的荧光粒子光斑几何光学特征相同，实际处理过程中极易出现错误跟踪，计算结果的准确性很难保证。


发明内容
针对上述问题，本发明的目的是提供一种具有较高的准确率、效率及高时间-空间分辨率的细胞原位主动变形测量方法。
为实现上述目的，本发明采取以下技术方案一种细胞原位主动变形测量方法，其包括以下步骤(1)分别在未发生弹性变形及发生弹性变形的细胞基底荧光图的参考图像F(x，y)和变形图像G(x，y)上，以设定的采样间距S为边长划分方形采样网格，单位为像素；(2)分别以各采样点(α，β)为中心，在参考图像F(x，y)上设定边长为(2N+1)的方形相关子区fα，β(x，y)，在变形图像G(x，y)上选择边长为(2M+1)的方形搜索区域Gα，β(x，y)，其中N的取值范围为15～30之间的整数，M的取值范围为20～80之间的整数，单位均为像素，且M＞N；(3)根据参考图像F(x，y)和变形图像G(x，y)的数字图像灰度分布结构表f(x，y)和g(x，y)，利用递推方法分别建立图像灰度全局求和结构表Sf、Sg及图像能量全局求和结构表


(4)根据数字图像灰度分布结构表f(x，y)和g(x，y)，计算以采样点(α，β)为中心的相关子区fα，β(x，y)与相应的搜索区域Gα，β(x，y)之间的互相关系数矩阵Pα，β，得到其快速递推关系式为Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)，其中，令

(u，v)表示参考图像的相关子区fα，β(x，y)中采样点(α，β)在变形图像的搜索区域Gα，β(x，y)中的位移；D是相关子区fα，β(x，y)的边长，D＝2N+1；(5)根据所述步骤(3)、(4)的结果，计算参考图像的相关子区fα，β(x，y)与相应变形图像的搜索区域Gα，β(x，y)的变形子区gα，β(x，y)之间的零均值归一化互相关系数矩阵Gα，β(u，v)；(6)采用基于梯度的亚像素位移定位算法，在零均值归一化互相关系数矩阵的最大峰值所在位置附近进行亚像素插值运算，求得参考图像中采样点(α，β)在变形图像中的准确位置为

式中，(Uα，β，Vα，β)表示采样点(α，β)在变形图像中的总位移大小，(uα，β，vα，β)是由所述零均值归一化相关系数矩阵Cα，β(u，v)中最大元素所确定的整像素位移，(Δuα，β，Δvα，β)是(uα，β，vα，β)相应的亚像素位移；(7)重复所述步骤(4)～步骤(6)，计算所有采样点在变形图像中的准确位置，进而得到整张变形图像位移场。
所述步骤(4)中，所述互相关系数矩阵Pα，β快速递推关系式由以下步骤得到①对于采样点(α，β)，令相关子区fα，β(x，y)与相应的搜索区域Gα，β(x，y)之间构成互相关系数矩阵为Pα，β；②与采样点(α，β)相邻的采样点(α+S，β)，其对应的相关子区fα+S，β(x，y)与相应的搜索区域Gα+S，β(x，y)之间的互相关系数矩阵为Pα+S，β；③将两个互相关系数矩阵Pα，β和Pα+S，β对比得到Pα+S，β(u，v)＝Pα，β(u，v)-Jα，β(u，v)+Iα，β(u，v)，其中，令

得到Jα，β(u，v)＝Iα-D，β(u，v)，D是相关子区fα，β(x，y)的边长，D＝2N+1；④由所述步骤①～③中各式能建立快速递推关系为Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)。
所述步骤(5)中，所述零均值归一化互相关系数矩阵Cα，β(u，v)为令其中，P(α，β；u，v)为图像参考子区fα，β(x，y)与相应搜索区域Gα，β(x，y)之间的互相关系数矩阵，其由所述步骤(4)中的所述互相关系数矩阵Pα，β快速递推关系式得到；F(α，β)为参考图像灰度平方和；G(α，β；u，v)为变形图像灰度平方和；Q(α，β；u，v)、F(α，β)和G(α，β；u，v)通过快速查询所述步骤(3)中预先建立的四个图像灰度、图像能量全局求和结构表得到。
本发明由于采取以上技术方案，其具有以下优点1、本发明由于采用了基于快速递推关系式和快速查询四个图像灰度、图像能量全局求和结构表的方式，快速计算零均值归一化相关系数矩阵，并根据该矩阵的峰值位置确定采样点的整像素位移大小，进而得到采样点的位移，因此实现了测量精确高效、高时间-空间分辨率的效果。2、本发明由于采用了快速查表方式和快速递推方法进行计算，因此可以有效地降低采样点位移计算的复杂性，特别是在高时间空间位移采样以及较大搜索区域尺寸情况下，能显著的提高细胞-基底位移场求解效率。3、本发明由于采用了快速递推方式进行计算，这种方法比传统数字图像相关方法的计算效率大约提高了10到50倍，特别是当采样网格边长减小时，计算效率尤为突出。因此，本发明可以应用于高时空分辨率的细胞原位主动变形测量中。



图1是本发明的整体流程示意图 图2是本发明采样点网格划分的相关子区、搜索区域示意图 图3是本发明的实施例一中单细胞粘附在聚丙烯酰胺弹性基底上的相差图 图4是本发明的实施例一中聚丙烯酰胺基底的变形前参考图像 图5是本发明的实施例一中聚丙烯酰胺基底的变形图像 图6是本发明的实施例一中基于本发明测量方法计算得到的荧光基底位移场
具体实施例方式 下面结合附图和实施例对本发明进行详细的描述。
如图1所示，本发明是基于零均值归一化互相关系数在细胞基底位移场求解过程中所表现出的内在特点，通过构建与图像灰度和图像能量密切相关的四个全局求和表，发展了一整套位移场求解快速递推方法，其步骤如下 1)将未发生弹性变形的细胞基底荧光图作为参考图像F(x，y)，发生弹性变形的细胞基底荧光图作为变形图像G(x，y)；设定采样间距为S，分别在参考图像F(x，y)和变形图像G(x，y)上以S为边长划分规则的方形采样网格(如图2所示)，单位为像素； 2)分别以每个采样点(α，β)为中心，在参考图像F(x，y)上设定边长为(2N+1)的方形相关子区fα，β(x，y)，在变形图像G(x，y)上选择边长为(2M+1)的方形搜索区域Gα，β(x，y)，其中N的取值范围一般为15～30之间的任意整数(其单位为像素)，而M的取值范围一般为20～80之间的任意整数(其单位也为像素)，且M＞N； 3)根据参考图像F(x，y)和变形图像G(x，y)的数字图像灰度分布结构表f(x，y)和g(x，y)(如表1、表2所示)，利用递推方法分别建立参考图像F(x，y)和变形图像G(x，y)的图像灰度全局求和结构表Sf、Sg及图像能量全局求和结构表


如表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 参考数字图像灰度分布结构表f(x，y) 表2 变形后数字图像灰度分布结构表g(x，y) 表3 参考图相灰度求和结构表Sf 表4 变形图相灰度求和结构表Sg 表5 参考图像能量求和结构表


表6 参考图像能量求和结构表

上述各表格中，p、q、w、H均为参考的图像参考子区fα，β(x，y)与相应搜索区域Gα，β(x，y)内任意点的取值，其单位为像素； 4)运用递推方法，根据参考图像F(x，y)和变形图像G(x，y)的数字图像灰度分布结构表f(x，y)和g(x，y)，计算以采样点(α，β)为中心的参考图像的相关子区fα，β(x，y)与相应变形图像的搜索区域Gα，β(x，y)之间的互相关系数矩阵Pα，β，进而得到互相关系数矩阵Pα，β快速递推关系式为 Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)，(5) 令 其中，(u，v)表示参考图像的相关子区fα，β(x，y)中采样点(α，β)在变形图像的搜索区域Gα，β(x，y)中的位移(如图2所示)；D是相关子区fα，β(x，y)的边长，D＝2N+1； 5)根据步骤3)和步骤4)的结果，计算参考图像的相关子区fα，β(x，y)与相应变形图像的搜索区域Gα，β(x，y)的变形子区gα，β(x，y)之间的零均值归一化互相关系数矩阵(Zero-normalized cross-correlation coefficient，ZNCC)，该零均值归一化互相关系数矩阵Cα，β(u，v)为 令 其中，P(α，β；u，v)为图像参考子区fα，β(x，y)与相应搜索区域Gα，β(x，y)之间的互相关系数矩阵；F(α，β)为参考图像灰度平方和；G(α，β；u，v)为变形图像灰度平方和； 公式(7)可由标准的零均值归一化互相关系数表达式经化简得到下式 上述公式中，

和

分别表示参考图像灰度和变形图像灰度的均值； 6)采用基于梯度的亚像素位移定位算法，在步骤5)中的ZNCC系数矩阵最大峰值所在位置附近进行亚像素插值运算，这样可以把荧光弹性基底的位移测量精度由整像素水平延伸至亚像素水平，即求得参考图像中采样点(α，β)在变形图像中的准确位置为 上式中，(Uα，β，Vα，β)表示采样点(α，β)在变形图像中的总位移大小，(uα，β，vα，β)是由步骤5)中的零均值归一化相关系数矩阵Cα，β(u，v)中最大元素所确定的整像素位移，(Δuα，β，Δvα，β)是(uα，β，vα，β)相应的亚像素位移； 7)重复步骤4)～步骤6)，计算所有采样点在变形后图像中的准确位置，进而得到整张变形图像位移场。
上述步骤4)中，互相关系数矩阵Pα，β快速递推关系式通过以下步骤得到 ①对于采样点(α，β)，令参考图像的相关子区fα，β(x，y)与相应的变形图像的搜索区域Gα，β(x，y)之间构成的互相关系数矩阵Pα，β为
②同理，对于与采样点(α，β)相邻的采样点(α+S，β)而言，其对应的相关子区fα+S，β(x，y)与相应的搜索区域Gα+S，β(x，y)之间的互相关系数矩阵Pα+S，β可以表示为
③将上述两个互相关系数矩阵Pα，β和Pα+S，β进行对比后，可以得到 其中，令 显然，将Jα，β(u，v)与公式(6)进行对比后可知，Jα，β(u，v)＝Iα-D，β(u，v)，D是相关子区fα，β(x，y)的边长，即D＝2N+1； ④由上述步骤中各式可建立如下快速递推关系 Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)。
上述步骤5)中，互相关系数矩阵P(α，β；u，v)根据步骤4)中的互相关系数矩阵Pα，β快速递推关系式可以得到 公式(7)中的Q(α，β；u，v)、F(α，β)和G(α，β；u，v)可以通过快速查询步骤3)中预先建立的四个图像灰度、图像能量全局求和结构表(表3～表6)得到，即 上述步骤6)中，亚像素位移(Δuα，β，Δvα，β)为 上式中，

和

分别表示在u和v方向上灰度的一阶梯度，灰度梯度算子取为[1/12，-8/12，0，8/12，-1/12]，这时 下面通过具体实施例对本发明的测量方法进行进一步的描述。
实施例一如图3～图6所示，给出的一组细胞变形图像，其图像尺寸为1027×533像素，当参考子区边长为40和60个像素时的计算效率测试结果如表7、表8所示， 表7 传统数字图像相关方法与本发明的测量方法计算时间之比(即加速比)
表8 传统数字图像相关方法与本发明的测量方法计算时间之比(即加速比)
从上述两个表中可以看到，当搜索半径介于10～50个像素，采样网格边长为5～20个像素时，传统数字图像相关方法与本发明的测量方法计算时间之比(即加速比)大约在10～50之间变化，这表明本发明的测量方法能够较传统数字图像相关方法提高计算效率大约10到50倍。特别是当采样网格边长减小时，计算效率尤为突出，这一特点预示着本发明的测量方法更适于高时空分辨率的细胞原位主动变形测量。
上述各实施例仅用于说明本发明，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。
权利要求
1.一种细胞原位主动变形测量方法，其包括以下步骤
(1)分别在未发生弹性变形及发生弹性变形的细胞基底荧光图的参考图像F(x，y)和变形图像G(x，y)上，以设定的采样间距S为边长划分方形采样网格，单位为像素；
(2)分别以各采样点(α，β)为中心，在参考图像F(x，y)上设定边长为(2N+1)的方形相关子区fα，β(x，y)，在变形图像G(x，y)上选择边长为(2M+1)的方形搜索区域Gα，β(x，y)，其中N的取值范围为15～30之间的整数，M的取值范围为20～80之间的整数，单位均为像素，且M＞N；
(3)根据参考图像F(x，y)和变形图像G(x，y)的数字图像灰度分布结构表f(x，y)和g(x，y)，利用递推方法分别建立图像灰度全局求和结构表Sf、Sg及图像能量全局求和结构表
(4)根据数字图像灰度分布结构表f(x，y)和g(x，y)，计算以采样点(α，β)为中心的参考图像的相关子区fα，β(x，y)与相应变形图像的搜索区域Gα，β(x，y)之间的互相关系数矩阵Pα，β，得到其快速递推关系式为
Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)，
其中，令
(u，v)表示参考图像的相关子区fα，β(x，y)中采样点(α，β)在变形图像的搜索区域Gα，β(x，y)中的位移；D是相关子区fα，β(x，y)的边长，D＝2N+1；
(5)根据所述步骤(3)、(4)的结果，计算参考图像的相关子区fα，β(x，y)与相应变形图像的搜索区域Gα，β(x，y)的变形子区gα，β(x，y)之间的零均值归一化互相关系数矩阵Cα，β(u，v)；
(6)采用基于梯度的亚像素位移定位算法，在零均值归一化互相关系数矩阵的最大峰值所在位置附近进行亚像素插值运算，求得参考图像中采样点(α，β)在变形图像中的准确位置为
式中，(Uα，β，Vα，β)表示采样点(α，β)在变形图像中的总位移大小，(uα，β，vα，β)是由所述零均值归一化相关系数矩阵Cα，β(u，v)中最大元素所确定的整像素位移，(Δuα，β，Δvα，β)是(uα，β，vα，β)相应的亚像素位移；
(7)重复所述步骤(4)～步骤(6)，计算所有采样点在变形图像中的准确位置，得到整张变形图像位移场。
2.如权利要求1所述的一种细胞原位主动变形测量方法，其特征在于所述步骤(4)中，所述互相关系数矩阵Pα，β快速递推关系式由以下步骤得到
①对于采样点(α，β)，令相关子区fα，β(x，y)与相应的搜索区域Gα，β(x，y)之间构成互相关系数矩阵为Pα，β；
②与采样点(α，β)相邻的采样点(α+S，β)，其对应的相关子区fα+S，β(x，y)与相应的搜索区域Gα+S，β(x，y)之间的互相关系数矩阵为Pα+S，β；
③将所述两个互相关系数矩阵Pα，β和Pα+S，β进行对比后，得到
Pα+S，β(u，v)＝Pα，β(u，v)-Jα，β(u，v)+Iα，β(u，v)
其中，令
D是相关子区fα，β(x，y)的边长，D＝2N+1；
④由所述步骤①～③中各式能建立快速递推关系为
Pα+S，β(u，v)＝Pα，β(u，v)-Iα-D，β(u，v)+Iα，β(u，v)。
3.如权利要求1所述的一种细胞原位主动变形测量方法，其特征在于所述步骤(5)中，所述零均值归一化互相关系数矩阵Cα，β(u，v)为
令
其中，P(α，β；u，v)为图像参考子区fα，β(x，y)与相应搜索区域Gα，β(x，y)之间的互相关系数矩阵，其由所述步骤(4)中的所述互相关系数矩阵Pα，β快速递推关系式得到；F(α，β)为参考图像灰度平方和；G(α，β；u，v)为变形图像灰度平方和；Q(α，β；u，v)、F(α，β)和G(α，β；u，v)通过快速查询所述步骤(3)中预先建立的四个图像灰度、图像能量全局求和结构表得到。
全文摘要
本发明涉及一种细胞原位主动变形测量方法，其测量步骤包括(1)分别在细胞-基底荧光参考图像和相应的变形图像上划分方形采样网格；(2)确定参考子区和搜索区域的尺寸；(3)建立四个全局求和表；(4)以递推方式计算以采样点为中心的图像参考子区与相应搜索区域之间的互相关系数；(5)利用步骤(4)的结果并通过查求和表的方式快速计算零均值归一化相关系数矩阵，并根据此矩阵的峰值位置确定采样点的整像素位移大小；(6)利用梯度法进行亚像素插值运算；(7)递推计算所有采样点处的位移，得到整张变形图像精确位移场。本发明具有测量精确高效、高时间-空间分辨率的效果。本发明可以应用于高时空分辨率的细胞原位主动变形测量中。
文档编号G01N33/48GK101813693SQ20101017069
公开日2010年8月25日 申请日期2010年5月6日 优先权日2010年5月6日
发明者黄建永, 潘晓畅, 秦雷, 朱涛, 熊春阳, 方竞 申请人:北京大学