一种显微图像序列中癌细胞轨迹追踪与关联方法与流程

文档序号:11865962阅读:316来源:国知局

本发明涉及细胞成像技术研究领域,尤其是一种癌细胞细胞轨迹跟踪与关联方法。



背景技术:

细胞轨迹追踪就是结合视频前后帧的信息以及细胞分裂生长的过程,实现对连续帧间细胞生命轨迹的追踪和关联。而本发明的目标对象是受药物抑制的癌细胞,其分裂生长过程较正常细胞具有特殊性,复杂性。传统的单一的细胞追踪方法并不适用。

从目前已有的多种目标跟踪方法来看,近年来研究较多的多目标的细胞跟踪主要分为两类。第一类基于帧间数据关联的算法,即通过前一帧信息来估计当前帧细胞情况,如此循环迭代。第二类是通过全局数据关联的多目标跟踪方法,利用更大的时间窗口获得更多的信息,通过全局轨迹长期的检测,在给定时间段内全局的解决轨迹连接问题。这些方法在处理树形层级轨迹关联时(即母细胞分裂到子细胞的情况)仍有一定的问题,需要我们在细胞分类后连续观测几帧细胞的情况以对细胞关联关系进行修正。



技术实现要素:

为了解决现有技术中一个目标对象的跟踪门限中量测数据可能不止一个(即量测数据可能来自于正确目标对象,也有可能来自于其他目标对象,也可能来自于杂波)的问题、精度较低,本发明提出了一种精度较高的显微图像序列中癌细胞轨迹追踪与关联方法。

本发明解决其技术问题所采用的技术方案是:

一种显微图像序列中癌细胞轨迹追踪与关联方法,所述细胞轨迹追踪与关联方法包括如下步骤:

1)基于最近邻法的癌细胞轨迹片段形成,即将检测到的癌细胞区域前后帧关联,形成可信的轨迹片段,过程如下:

首先,设置跟踪波门限r,落到跟踪波门中的量测作为候选回波,即目标的量测值zi(k)是否能够满足公式:

<mrow> <msup> <mrow> <mo>&lsqb;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>z</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mi>T</mi> </msup> <msup> <mi>S</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&lsqb;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>z</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> <mo>&le;</mo> <mi>r</mi> </mrow>

其中是跟踪波门的中心,如果波门中只有一个量测,该量测值则直接被用于轨迹的更新中,形成连续的轨迹段;如若有一个以上的候选回波,则应根据距离计算方式,找到距离最近的候选回波,来用于轨迹的更新;其中,zi(k)对应的统计距离为:

<mrow> <msup> <mi>d</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>(</mo> <mi>k</mi> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>&lsqb;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>z</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mi>T</mi> </msup> <msup> <mi>S</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&lsqb;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>z</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow>

接着,考虑当跟踪波门中只有一个量测,或是距离最近的候选回波距离前一帧目标的距离d足够近的情况下,将其连接到细胞轨迹中,形成细胞轨迹片段;而当最近量测值距离超过某一阈值,或是有两个及以上候选回波距离前一帧目标的距离d都比较近的情况下,细胞轨迹片段断裂,重新开始新的轨迹;

2)基于细胞动态特征的全局轨迹关联;

由于细胞在生长、分裂、融合的过程中形态不断发生改变,仅依靠单一的特征难以正确区分,因此结合细胞的动态特征匹配,将再次断裂的细胞轨迹片段关联。

细胞动态特征参数包括以下4种:运动位移参数Edisplacement、运动偏斜参数Eskewness、面积变化参数Earea以及形状变化参数Edeformation

细胞u和v的运动位移参数Edisplacement定义如下:

<mrow> <msub> <mi>E</mi> <mrow> <mi>d</mi> <mi>i</mi> <mi>s</mi> <mi>p</mi> <mi>l</mi> <mi>a</mi> <mi>c</mi> <mi>e</mi> <mi>m</mi> <mi>e</mi> <mi>n</mi> <mi>t</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <msubsup> <mi>c</mi> <mi>u</mi> <mi>n</mi> </msubsup> <msubsup> <mi>c</mi> <mi>v</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>|</mo> <mo>|</mo> </mrow> <msqrt> <mrow> <msup> <mi>H</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>W</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow>

其中,cu,cv分别是细胞u、v在第n帧和第n+1帧的位置,H和W分别指当前帧图像的长和宽;

运动偏斜参数Eskewness定义如下:

这里ct,cu,cv分别是细胞t,u,v在第n-1帧,第n帧以及第n+1帧时的中心位置;这个参数是用了测量运动方向的偏移。

面积变化参数Earea用于测量细胞u和v的重叠程度,定义如下:

<mrow> <msub> <mi>E</mi> <mrow> <mi>a</mi> <mi>r</mi> <mi>e</mi> <mi>a</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>{</mo> <mn>1</mn> <mo>-</mo> <mfrac> <msubsup> <mi>s</mi> <mi>c</mi> <mrow> <mi>n</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <msup> <mn>1</mn> <mn>2</mn> </msup> </mrow> </msubsup> <mrow> <msubsup> <mi>s</mi> <mi>u</mi> <mi>n</mi> </msubsup> <msubsup> <mi>s</mi> <mi>v</mi> <mrow> <mi>n</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> </mrow> </mfrac> <mo>}</mo> </mrow>

其中,指细胞u在第n帧的面积,值细胞v在第n+1帧中的面积,指两细胞的重叠程度;这个公式表明,两细胞重叠面积越多,两者就更倾向于是同一个细胞。

对于形状变化参数参数Edeformation,使用匹配细胞u和v面积的椭圆离心率,椭圆离心率公式Q=P2/(4π×A2),这里P和A分别表示椭圆的周长和面积,形变参数Edeformation定义如下:

<mrow> <msub> <mi>E</mi> <mrow> <mi>d</mi> <mi>e</mi> <mi>f</mi> <mi>o</mi> <mi>r</mi> <mi>m</mi> <mi>a</mi> <mi>t</mi> <mi>i</mi> <mi>o</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mrow> <msubsup> <mi>Q</mi> <mi>u</mi> <mi>n</mi> </msubsup> <mo>-</mo> <msubsup> <mi>Q</mi> <mi>v</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> </mrow> <mo>|</mo> </mrow> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msubsup> <mi>Q</mi> <mi>u</mi> <mi>n</mi> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msubsup> <mi>Q</mi> <mi>u</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow>

使用4个参数来估计轨迹片段间的相似程度,另外,使用一个基于双边权值匹配的全局最优技术来最小化总代价,函数定义如下:

E(ci,cj)=α1Edisplacement2Eskewness3Earea4Edeformation

其中,E(ci,cj)就是第n帧中的细胞i与第n+1帧中的细胞j的代价函数,αi限定在0到1之间,且

通过计算各个细胞的动态特征,得到代价函数的值之后,值最小的就是正确匹配的细胞,即可将此处断裂的细胞轨迹片段正确关联起来。

3)基于分类检测的轨迹关联,过程如下:

在数据关联中,无论是基于最近邻法的癌细胞轨迹片段形成还是基于动态特征的全局轨迹关联,都是一对一的匹配过程。如若涉及到细胞分裂的情况,则需要根据分裂的具体状况进行分类检测以关联轨迹。

首先,根据细胞分裂时候的特征判断是否属于分裂事件;

再者,因为本发明研究对象为受药物抑制的癌细胞,正常细胞的分裂检测事件判断对本数据集并不适用。根据癌细胞分裂特征,判断是否将未关联的轨迹片段关联到母细胞,并标记分裂事件。

4)基于轨迹关联的状态检测纠正;

在癌细胞生命轨迹关联之后,即观测细胞都已形成连续的轨迹后,根据前后帧细胞的状态对于误检测的细胞状态进行修正;

从第一帧开始,记录下待监测癌细胞所处的阶段,在该细胞的整个轨迹过程中,向后逐步检测,记录下与前一帧状态不同的时刻;若所记录的时刻为突变点,即相对整个序列的异常点,则根据前后细胞状态修正误检测的情况;

另外,由于细胞生长遵循分裂间期、分裂中期和分裂后期这四个阶段,细胞轨迹中的状态一定也是按此顺序出现,由于各个癌细胞所处阶段不同,起始状态可能是分裂间期、分裂前期或是分裂中期,但在分裂后期之前,都遵循正常的规律。分裂后期之后,即受药物抑制癌细胞挣扎后的状态可能会有所不同;

遵循如上两个原则以修正癌细胞细胞跟踪轨迹。

本发明的有益效果主要表现在:基于数据关联的癌细胞细胞轨迹追踪方法,既能实现对癌细胞细胞轨迹的关联,同时利用分类检测和轨迹关联的状态检测又能更好的关联轨迹片段及对其加以修正。

具体实施方式

下面对本发明作进一步描述。

一种显微图像序列中癌细胞轨迹追踪与关联方法,所述追踪与关联方法包括如下步骤:

1)基于最近邻法的癌细胞轨迹片段形成,即将检测到的癌细胞区域前后帧关联,形成可信的轨迹片段,过程如下:

首先,设置跟踪波门限r,落到跟踪波门中的量测作为候选回波,即目标的量测值zi(k)是否能够满足公式:

<mrow> <msup> <mrow> <mo>&lsqb;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>z</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mi>T</mi> </msup> <msup> <mi>S</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&lsqb;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>z</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> <mo>&le;</mo> <mi>r</mi> </mrow>

其中是跟踪波门的中心,如果波门中只有一个量测,该量测值则直接被用于轨迹的更新中,形成连续的轨迹段;如若有一个以上的候选回波,则应根据距离计算方式,找到距离最近的候选回波,来用于轨迹的更新;其中,zi(k)对应的统计距离为:

<mrow> <msup> <mi>d</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>(</mo> <mi>k</mi> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>&lsqb;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>z</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mi>T</mi> </msup> <msup> <mi>S</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&lsqb;</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <mi>z</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow>

接着,考虑当跟踪波门中只有一个量测,或是距离最近的候选回波距离前一帧目标的距离d足够近的情况下,将其连接到细胞轨迹中,形成细胞轨迹片段。而当最近量测值距离超过某一阈值,或是有两个及以上候选回波距离前一帧目标的距离d都比较近的情况下,细胞轨迹片段断裂,重新开始新的轨迹;

2)基于细胞动态特征的全局轨迹关联;

由于细胞在生长、分裂、融合的过程中形态不断发生改变,仅依靠单一的特征难以正确区分,因此结合细胞的动态特征匹配,将再次断裂的细胞轨迹片段关联;

细胞动态特征参数包括以下4种:运动位移参数(Edisplacement)、运动偏斜参数(Eskewness)、面积变化参数(Earea)、以及形状变化参数(Edeformation)。

细胞u和v的运动位移参数Edisplacement定义如下:

<mrow> <msub> <mi>E</mi> <mrow> <mi>d</mi> <mi>i</mi> <mi>s</mi> <mi>p</mi> <mi>l</mi> <mi>a</mi> <mi>c</mi> <mi>e</mi> <mi>m</mi> <mi>e</mi> <mi>n</mi> <mi>t</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <msubsup> <mi>c</mi> <mi>u</mi> <mi>n</mi> </msubsup> <msubsup> <mi>c</mi> <mi>v</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>|</mo> <mo>|</mo> </mrow> <msqrt> <mrow> <msup> <mi>H</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>W</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow>

其中,cu,cv分别是细胞u、v在第n帧和第n+1帧的位置,H和W分别指当前帧图像的长和宽。

运动偏斜参数Eskewness定义如下:

其中,ct,cu,cv分别是细胞t,u,v在第n-1帧,第n帧以及第n+1帧时的中心位置,运动偏斜参数是用了测量运动方向的偏移;

面积变化参数Earea用于测量细胞u和v的重叠程度,定义如下:

<mrow> <msub> <mi>E</mi> <mrow> <mi>a</mi> <mi>r</mi> <mi>e</mi> <mi>a</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>{</mo> <mn>1</mn> <mo>-</mo> <mfrac> <msubsup> <mi>s</mi> <mi>c</mi> <mrow> <mi>n</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <msup> <mn>1</mn> <mn>2</mn> </msup> </mrow> </msubsup> <mrow> <msubsup> <mi>s</mi> <mi>u</mi> <mi>n</mi> </msubsup> <msubsup> <mi>s</mi> <mi>v</mi> <mrow> <mi>n</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> </mrow> </mfrac> <mo>}</mo> </mrow>

其中,指细胞u在第n帧的面积,值细胞v在第n+1帧中的面积,指两细胞的重叠程度。这个公式表明,两细胞重叠面积越多,两者就更倾向于是同一个细胞。

对于形状变化参数Edeformation,使用匹配细胞u和v面积的椭圆离心率,椭圆离心率公式Q=P2/(4π×A2),这里P和A分别表示椭圆的周长和面积,形状变化参数Edeformation定义如下:

<mrow> <msub> <mi>E</mi> <mrow> <mi>d</mi> <mi>e</mi> <mi>f</mi> <mi>o</mi> <mi>r</mi> <mi>m</mi> <mi>a</mi> <mi>t</mi> <mi>i</mi> <mi>o</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mrow> <msubsup> <mi>Q</mi> <mi>u</mi> <mi>n</mi> </msubsup> <mo>-</mo> <msubsup> <mi>Q</mi> <mi>v</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> </mrow> <mo>|</mo> </mrow> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msubsup> <mi>Q</mi> <mi>u</mi> <mi>n</mi> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msubsup> <mi>Q</mi> <mi>u</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow>

使用4个参数来估计轨迹片段间的相似程度,另外,使用一个基于双边权值匹配的全局最优技术来最小化总代价,函数定义如下:

E(ci,cj)=α1Edisplacement2Eskewness3Earea4Edeformation

这里的E(ci,cj)就是第n帧中的细胞i与第n+1帧中的细胞j的代价函数。这里αi限定在0到1之间,且

通过计算各个细胞的动态特征,得到代价函数的值之后,值最小的就是正确匹配的细胞,即可将此处断裂的细胞轨迹片段正确关联起来;

3)基于分类检测的轨迹关联

在数据关联中,无论是基于最近邻法的癌细胞轨迹片段形成还是基于动态特征的全局轨迹关联,都是一对一的匹配过程。如若涉及到细胞分裂的情况,则需要根据分裂的具体状况进行分类检测以关联轨迹。

首先,根据细胞分裂时候的特征判断是否属于分裂事件,特征如下:

正常细胞发生分裂时,在下一刻会有两个独立的细胞产生,这两个细胞距离前一时刻的母细胞距离较近;

母细胞在一份为二之前,细胞形状变圆,周围发亮;

两个子细胞之间相距很近;

独立的两个子细胞面积相对较小,二者面积和近似于前一帧中母细胞的面积。

再者,因为本发明研究对象为受药物抑制的癌细胞,正常细胞的分裂检测事件判断对本数据集并不适用。首先,癌细胞若发生分裂,未必是一分为二,可能分裂为两个或是两个以上子细胞;其次,母细胞形状变圆,变圆发亮这个阶段,药物抑制下的癌细胞也会经历(即第三个阶段),但通常在这一阶段后癌细胞会与药物对抗做挣扎,若有分裂也是在挣扎失败之后发生,因此变圆发亮这一特征无法作为判别依据;此外,癌细胞在受药物作用挣扎的过程中运动剧烈,与之前相比可能有较大的位移偏差。根据癌细胞的这些分裂特征,对于在时空序列中途出现的未能与前帧关联起来的癌细胞轨迹片段在前一帧中该细胞附近寻找是否有处于第四阶段的细胞存在,若有,将未关联的轨迹片段关联到母细胞,并标记分裂事件。若未找到合适的匹配,则认为该细胞为由外部进入的新出现细胞;

4)基于轨迹关联的状态检测纠正

在癌细胞生命轨迹关联之后,即观测细胞都已形成连续的轨迹后,根据前后帧细胞的状态对于误检测的细胞状态进行修正。

从第一帧开始,记录下待监测癌细胞所处的阶段,在该细胞的整个轨迹过程中,向后逐步检测,记录下与前一帧状态不同的时刻。若所记录的时刻为突变点,即相对整个序列的异常点,则根据前后细胞状态修正误检测的情况。

另外,由于细胞生长遵循分裂间期、分裂中期和分裂后期这四个阶段,细胞轨迹中的状态一定也是按此顺序出现。由于各个癌细胞所处阶段不同,起始状态可能是分裂间期、分裂前期或是分裂中期,但在分裂后期之前,都遵循正常的规律。分裂后期之后,即受药物抑制癌细胞挣扎后的状态可能会有所不同。

遵循如上两个原则,选取异常检测方法,如常见的基于统计的方法,基于距离的方法,基于密度的方法,基于偏离的方法,基于聚类的方法等以修正癌细胞细胞跟踪轨迹。

当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1