一种运动粒子的二维及三维示踪轨迹续连方法与流程

1.本发明涉及微观和宏观粒子示踪技术领域，特别涉及一种运动粒子的二维及三维示踪轨迹续连方法。


背景技术：

2.粒子示踪技术在流体力学、气象学、化学、生物医学、微生物学、海洋科学等领域均有广泛应用。它可用于研究流场的三维速度分布、宏观生物动植物的运动行为、各类细胞、微生物和颗粒的微观行为和动力学过程、海洋污染物的示踪溯源等研究课题，是搭载于商用光学显微镜等成像仪器上的常见技术。粒子示踪技术可分为二维示踪和三维示踪，所示踪的目标尺寸从宏观到纳米级不等。目前的粒子定位技术普遍存在鲁棒性差，易受到环境噪声干扰的缺点，导致示踪中断，轨迹错判或定位偏差等系统性误差的出现，严重限制了粒子示踪技术的发展和广泛应用。
3.粒子示踪技术是一种高度专门化的技术，在不同研究领域，需要使用不同的定位、示踪、优化方法。定位方法主要包括二维图像识别、双目视觉、全息技术等，示踪方法主要包括最近邻法、位置预测法、贪婪算法、全局优化法等。这些方法都不能完全描述和处理粒子的复杂运动，对于三维环境和高浓度粒子的示踪效果则相较二维示踪方法更差，轨迹可发生频繁断连，导致轨迹碎片化，无法对示踪目标的运动状态得出统计意义上的描述。因此，建立一种用于修复粒子示踪过程中因为环境变化和粒子相互干扰影响导致的轨迹中断，恢复丢失信息的轨迹续连技术是极有必要的。该技术具有普适性，可解决轨迹断连这一粒子示踪领域的基础问题，在生命科学、医学、流体力学、环境监测、工业生产等领域将发挥重要的作用。willneff等人开发并实现了一种新的时空匹配方法，并在不同的数据集上进行了测试。在正确利用冗余信息的情况下，该算法可以实现轨迹的连续。该方法实现了高通量的三维粒子跟踪，但定位精度仅为微米级，且粒子浓度增加时定位精度会降低。patel等人(patel,m.；leggett,s.e.；landauer,a.k.；wong,i.y.；franck,c.,rapid,topology-based particle tracking for high-resolution measurements of large complex 3d motion fields.sci rep 2018,8(1),5581.)提出了一种新的粒子跟踪方案，称为基于拓扑的粒子跟踪(t-pt)。这种方法可以高效率而且准确地在大空间和高空间梯度下跟踪大量粒子，但存在正确率较低，难以跟踪运动步长较大的颗粒的问题。


技术实现要素：

4.本发明的主要目的在于克服现有技术的缺点与不足，提供一种运动粒子的二维及三维示踪轨迹续连方法。
5.本发明至少通过如下技术方案之一实现。
6.一种运动粒子的二维及三维示踪轨迹续连方法，包括以下步骤：
7.步骤s1、选择前轨迹和后轨迹；
8.步骤s2、线性拟合前轨迹和后轨迹之间的丢失区域；
9.步骤s3、根据需求向前线性拟合前轨迹，向后线性拟合后轨迹到所需长度；
10.步骤s4、通过图像识别算法获得拟合位置样品图像中的粒子二维坐标；
11.步骤s5、通过位置预测算法筛选粒子二维或三维坐标；
12.步骤s6、若为三维轨迹则恢复拟合位置的三维信息；
13.步骤s7、验证前轨迹和后轨迹是否属于同一个粒子。
14.优选地，前轨迹l1为随机选择，后轨迹l2为前轨迹l1周围同时满足时间和空间双重限制条件的轨迹；
15.时间限制条件为1≤f
21-f
1end
≤f，f
1end
到f
21
之间的帧称为丢失区域，f为时间限制参数，f
1end
为前轨迹l1的最后一帧，f
21
为后轨迹l2的第一帧；
16.空间限制条件为0≤d
gap
≤maxd0 (f
21-f
1end-1)maxd1；其中，maxd0为单帧最大位移maxd0，maxd1为多帧最大位移，前轨迹l1的最后一个点到后轨迹l2的第一个点的距离为d
gap
；
17.对于任意的前轨迹l1，符合时间和空间限制条件的所有轨迹都能够选为属于该前轨迹l1的后轨迹l2。
18.优选地，步骤s2包括：将前轨迹l1和属于这条前轨迹l1的一条后轨迹l2进行线性拟合，用线段l连接前轨迹l1的最后一个点和后轨迹l2的第一个点，在这条线段l上均匀地取f
21-f
1end
1个位置，并去掉线段l两个端点的位置，即为前轨迹l1和后轨迹l2之间丢失区域的线性拟合位置，f
21
为后轨迹l2的第一帧，f
1end
为前轨迹l1的最后一帧。
19.优选地，步骤s3包括以下步骤：
20.设n为延长拟合参照帧数，n≥1并且小于前轨迹l1和后轨迹l2包含的帧数，设延长帧数m；
21.取前轨迹l1的前n个点，对这n个点的时间坐标t和坐标(x、y)进行线性回归，获得坐标(x、y)相对于时间坐标t的线性方程：
22.x(t)＝k
x
t b
x
23.y(t)＝kyt by24.其中k
x
为x相对于t的斜率，b
x
为x相对于t的截距，ky为y相对于t的斜率，by为y相对于t的截距；
25.设前轨迹l1的第一帧为f
11
，将线性方程x(t)、y(t)代入时间坐标t＝f
11-1到f
11-m即得前轨迹l1向前一帧到m帧的线性预测位置的坐标(x、y)；
26.若为三维轨迹，对这n个点的时间坐标t和坐标z进行线性回归，获得坐标z相对于时间坐标t的线性方程：
27.z(t)＝kzt bz28.其中kz为z相对于t的斜率，bz为z相对于t的截距；
29.设前轨迹l1的第一帧为f
11
，将线性方程z(t)代入时间坐标t＝f
11-1到f
11-m即得前轨迹l1向前一帧到m帧的线性预测位置的坐标z；
30.设n为延长拟合参照帧数，n≥1并且小于前轨迹l1和后轨迹l2包含的帧数，设延长帧数m；
31.取后轨迹l2的后n个点，对这n个点的时间坐标t和坐标(x、y)进行线性回归，获得坐标(x、y)相对于时间坐标t的线性方程：
32.x(t)＝k
x
t b
x
33.y(t)＝kyt by34.其中k
x
为x相对于t的斜率，b
x
为x相对于t的截距，ky为y相对于t的斜率，by为y相对于t的截距；
35.设后轨迹l2的最后一帧为f
2end
，将线性方程x(t)、y(t)代入时间坐标t＝f
2end
1到f
2end
m即得后轨迹l2向后一帧到m帧的线性预测位置的坐标(x、y)；
36.若为三维轨迹，对这n个点的时间坐标t和坐标z进行线性回归，获得坐标z相对于时间坐标t的线性方程：
37.z(t)＝kzt bz38.其中kz为z相对于t的斜率，bz为z相对于t的截距；
39.设后轨迹l2的最后一帧为f
2end
，将线性方程z(t)代入时间坐标t＝f
2end
1到f
2end
m即得后轨迹l2向后一帧到m帧的线性预测位置的坐标z；
40.步骤s2和步骤s3中获得的线性拟合位置统称为拟合位置，记为(xf，yf)，若为三维轨迹则记为(xf，yf，zf)。
41.优选地，所述图像识别算法包括以下内容：
42.设光强阈值为t，在样品图像中，选出所有光强大于t的像素，其中所有相邻的像素属于同一个粒子，将粒子包含的所有像素的二维坐标平均值或光强加权后的坐标平均值作为粒子的二维坐标(x
′
，y
′
)。
43.优选地，所述位置预测算法为：
44.设j为位置预测参照帧数，取前轨迹l1的后j个点，对这j个点的时间坐标t和坐标(x、y)进行线性回归，获得坐标(x、y)相对于时间坐标t的线性方程：
45.x(t)＝k
x
t b
x
46.y(t)＝kyt by47.其中k
x
为x相对于t的斜率，b
x
为x相对于t的截距，ky为y相对于t的斜率，by为y相对于t的截距；
48.将线性方程x(t)、y(t)代入t＝f
1end
1即得前轨迹l1向后一帧的线性预测位置的坐标(x、y)，获得的预测坐标记为(x
p
，y
p
)，其中f
1end
为前轨迹l1的最后一帧；
49.在前轨迹l1向后一帧中，计算所有粒子的二维坐标(x
′
，y
′
)到预测坐标(x
p
，y
p
)的距离，距离(x
p
，y
p
)最近的(x
′
，y
′
)选为前轨迹l1向后一帧的粒子二维坐标，记为(xr，yr)，把(xr，yr)置于前轨迹l1的末尾；
50.若为三维轨迹则把(xr，yr，zf)置于前轨迹l1的末尾，成为新的前轨迹l1的最后一点，重复该步骤，直到时间坐标t为丢失区域的中间帧；
51.设j为位置预测参照帧数，取后轨迹l2的前j个点，对这j个点的时间坐标t和坐标x、y进行线性回归，获得坐标x、y相对于时间坐标t的线性方程：
52.x(t)＝k
x
t b
x
53.y(t)＝kyt by54.其中k
x
为x相对于t的斜率，b
x
为x相对于t的截距，ky为y相对于t的斜率，by为y相对于t的截距；
55.将线性方程x(t)、y(t)代入t＝f
21-1即得后轨迹l2向前一帧的线性预测位置的坐标x、y，其中f
21
为后轨迹l2的第一帧，这样获得的预测坐标记为(x
p
，y
p
)；
56.在后轨迹l2向前一帧中，计算所有粒子的二维坐标(x
′
，y
′
)到预测坐标(x
p
，y
p
)的距离，距离(x
p
，y
p
)最近的(x
′
，y
′
)选为后轨迹l2向前一帧的粒子二维坐标，记为(xr，yr)，把(xr，yr)置于后轨迹l2的开头；
57.若为三维轨迹则把(xr，yr，zf)置于后轨迹l2的开头，成为新的后轨迹l2的第一点，重复该步骤，直到时间坐标t为丢失区域的中间帧的后一帧；
58.对前轨迹l1的向前线性拟合帧和后轨迹l2的向后线性拟合帧同样如此执行该步骤。
59.优选地，步骤s6具体包括：
60.若样品图像中含有三维信息，则在样品图像的坐标(xr，yr)处，使用该三维信息对应的方式获取(xr，yr，zf)的z坐标，记为zr，则该粒子坐标为(xr，yr，zr)；
61.对光场中每一层不同的z，其中坐标(xr，yr)对应的粒子所包含的所有像素的强度的和为iz，每个不同的z都有一个iz，对iz和z进行高斯拟合，拟合结果的最大值记为zr，粒子的三维坐标即为(xr，yr，zr)。
62.优选地，样品图像中含有三维信息为全息图像、双目视觉、全内反射光强图。
63.优选地，当三维信息为全息图像时，通过高斯拟合获得粒子的三维坐标：
64.设zg为高斯拟合范围，在坐标(xr，yr，zf)处，通过rayleigh-sommerfeld算法，对全息图像进行三维重建，获取z
f-zg到zf zg的三维光场，计算公式如下：
[0065][0066]
u(r，z)＝ft-1
(ft(is(r，0)
·
h(q，-z)))
[0067]
式中，h(r，-z)为瑞利索末菲传播子，r和z分别为拟合位置的横向坐标和纵向坐标，i为虚数单位，k为波数，r为光传播的距离，u(r，z)是三维重建所得的光场，is(r，0)为全息图像在r处的光强，ft-1
、ft分别为傅里叶逆变换和傅里叶变换，h(q，-z)为h(r，-z)的傅里叶变换，q为频域坐标。
[0068]
优选地，步骤s7的验证为：
[0069]
设纠错距离为wrong_d，计算前轨迹l1的最后一点到后轨迹l2的第一点的距离d
finish
，若大于wrong_d，则前轨迹l1和后轨迹l2不属于同一个粒子，否则前轨迹l1和后轨迹l2属于同一个粒子。
[0070]
本发明同时适用于二维或三维轨迹，能续连已有的碎片式轨迹，补齐碎片轨迹之间的丢失区域，若轨迹长度不足则向两端延长轨迹，有效增加了轨迹点数目。因此，该发明可显著提升粒子示踪技术的准确性和实用性，适用于对各种宏观和微观运动中的多个目标粒子进行同时示踪，包括但不限于细菌、细胞、流场中的微纳气泡和液滴、海面漂浮物、微小动植物等。
[0071]
与现有的技术相比，本发明的有益效果为：
[0072]
能续连已有的碎片式轨迹，补齐碎片轨迹之间的丢失区域，若轨迹长度不足则向两端延长轨迹，有效增加了轨迹点数目，增加了所获得的信息量，显著提升了粒子示踪技术的效率，准确性和实用性。
附图说明
[0073]
图1是本发明所述一种运动粒子的二维及三维示踪轨迹续连方法的方法流程图；
[0074]
图2是本发明实施例中选取后轨迹l2的时间限制示意图；
[0075]
图3是本发明实施例中选取后轨迹l2的空间限制示意图；
[0076]
图4是本发明实施例中线性拟合丢失区域的示意图；
[0077]
图5是本发明实施例中向两端线性拟合l1和l2的示意图；
[0078]
图6是本实施例中通过图像识别算法获得样品图像中的粒子二维坐标的示意图；
[0079]
图7是本发明实施例中通过位置预测算法筛选粒子二维坐标的示意图；
[0080]
图8是本发明实施例中粒子强度求和范围的示意图；
[0081]
图9是本发明实施例中通过高斯拟合获得粒子的三维坐标的示意图；
[0082]
图10是本发明实施例中验证前轨迹l1和后轨迹l2是否属于同一个粒子的示意图；
[0083]
图11是本发明实施例中两条碎片轨迹续连前后的示意图。
具体实施方式
[0084]
下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
[0085]
一种运动粒子的二维及三维示踪轨迹续连方法，如图1所示，包括以下步骤：
[0086]
第一步、选取前轨迹l1和后轨迹l2。所述轨迹包括常用的光学显微镜和摄像技术、双目视觉、数字全息显微镜、全内反射显微镜等各类粒子示踪技术获得的粒子二维或三维轨迹。
[0087]
前轨迹l1为随机选择，每条轨迹都会被选为前轨迹l1，且只会被选为前轨迹l1一次。后轨迹l2为前轨迹l1附近同时满足时间和空间双重限制条件的轨迹。
[0088]
时间限制条件为1≤f
21-f
1end
≤f，f
1end
到f
21
之间的帧称为丢失区域，f为时间限制参数，f
1end
为前轨迹l1的最后一帧，f
21
为后轨迹l2的第一帧，如图2所示，f
1end
到f
21
之间的帧(不含两端)称为丢失区域。
[0089]
如图3所示，空间限制条件为0≤d
gap
≤maxd0 (f
21-f
1end-1)maxd1，maxd0为单帧最大位移，maxd1为多帧最大位移，d
gap
为前轨迹l1的最后一个点到后轨迹l2的第一个点的距离对于一条前轨迹l1，符合时间和空间限制条件的所有轨迹都会选为属于该前轨迹l1的后轨迹l2。
[0090]
第二步、线性拟合l1和后轨迹l2之间的丢失区域，通过步骤s1获取的一条前轨迹l1可能有多条后轨迹l2。对一条前轨迹l1和属于这条前轨迹l1的一条后轨迹l2进行线性拟合，用一条线段连接前轨迹l1的最后一个点和后轨迹l2的第一个点。在这条线段上均匀地取f
21-f
1end
1个位置，并去掉端点的两个位置，即为前轨迹l1和后轨迹l2之间丢失区域的线性拟合位置。用此方法获取所有前轨迹l1和后轨迹l2的丢失区域的线性拟合位置，如图4所示。
[0091]
第三步、根据需求向前线性拟合前轨迹l1，向后线性拟合后轨迹l2到所需长度。
[0092]
设n为延长拟合参照帧数，n≥1并且小于前轨迹l1和后轨迹l2包含帧数。设延长帧数m。
[0093]
取前轨迹l1的前n个点，对这n个点的时间坐标t和坐标x、y进行线性回归，获得坐标x、y相对于时间坐标t的线性方程：
[0094]
x(t)＝k
x
t b
x
[0095]
y(t)＝kyt by[0096]
其中k
x
为x相对于t的斜率，b
x
为x相对于t的截距，ky为y相对于t的斜率，by为y相对于t的截距。
[0097]
设前轨迹l1的第一帧为f
11
，将线性方程x(t)、y(t)代入时间坐标t＝f
11-1到f
11-m即得前轨迹l1向前一帧到m帧的线性预测位置的坐标x、y，如图5所示；
[0098]
若为三维轨迹，对这n个点的时间坐标t和坐标z进行线性回归，获得坐标z相对于时间坐标t的线性方程：
[0099]
z(t)＝kzt bz[0100]
其中kz为z相对于t的斜率，bz为z相对于t的截距。
[0101]
设前轨迹l1的第一帧为f
11
，将线性方程z(t)代入时间坐标t＝f
11-1到f
11-m即得前轨迹l1向前一帧到m帧的线性预测位置的坐标z；
[0102]
设n为延长拟合参照帧数，n≥1并且小于前轨迹l1和后轨迹l2包含的帧数，设延长帧数m；延长帧数m可根据需求任意设置，但不宜过大，以免粒子实际运动轨迹偏离线性方向而无法进行图像识别。
[0103]
取后轨迹l2的后n个点，对这n个点的时间坐标t和坐标x、y进行线性回归，获得坐标x、y相对于时间坐标t的线性方程：
[0104]
x(t)＝k
x
t b
x
[0105]
y(t)＝kyt by[0106]
其中k
x
为x相对于t的斜率，b
x
为x相对于t的截距，ky为y相对于t的斜率，by为y相对于t的截距。
[0107]
设后轨迹l2的最后一帧为f
2end
，将线性方程x(t)、y(t)代入时间坐标t＝f
2end
1到f
2end
m即得后轨迹l2向后一帧到m帧的线性预测位置的坐标x、y；
[0108]
若为三维轨迹，对这n个点的时间坐标t和坐标z进行线性回归，获得坐标z相对于时间坐标t的线性方程：
[0109]
z(t)＝kzt bz[0110]
其中kz为z相对于t的斜率，bz为z相对于t的截距。
[0111]
设后轨迹l2的最后一帧为f
2end
，将线性方程z(t)代入时间坐标t＝f
2end
1到f
2end
m即得后轨迹l2向后一帧到m帧的线性预测位置的坐标z；
[0112]
步骤s2和步骤s3中获得的线性拟合位置统称为拟合位置，记为(xf，yf)。若为三维轨迹则记为(xf，yf，zf)。
[0113]
第四步、通过图像识别算法获得样品图像中的粒子二维坐标。
[0114]
设光强阈值为t，在样品图像中，选出所有光强大于t的像素。其中所有相邻的像素属于同一个粒子。一个粒子的二维坐标(x
′
，y
′
)是这个粒子包含的所有像素的二维坐标平均值或光强加权后的坐标平均值，如图6所示。由于样品图像中可能包含多个粒子，以及噪音的干扰，一个样品图像中可能获得多个(x
′
，y
′
)。
[0115]
第五步、通过位置预测算法筛选粒子二维坐标。
[0116]
设j为位置预测参照帧数，取前轨迹l1的后j个点，对这j个点的时间坐标t和坐标x、y进行线性回归，获得坐标x、y相对于时间坐标t的线性方程：
[0117]
x(t)＝k
x
t b
x
[0118]
y(t)＝kyt by[0119]
其中k
x
为x相对于t的斜率，b
x
为x相对于t的截距，ky为y相对于t的斜率，by为y相对于t的截距。
[0120]
将线性方程x(t)、y(t)代入t＝f
1end
1即得前轨迹l1向后一帧的线性预测位置的坐标x、y，其中f
1end
为前轨迹l1的最后一帧，这样获得的预测坐标记为(x
p
，y
p
)。
[0121]
在前轨迹l1向后一帧中，计算所有粒子的二维坐标(x
′
，y
′
)到预测坐标(x
p
，y
p
)的距离，距离(x
p
，y
p
)最近的(x
′
，y
′
)选为前轨迹l1向后一帧的粒子二维坐标，记为(xr，yr)。把(xr，yr)置于前轨迹l1的末尾，若为三维轨迹则把(xr，yr，zf)置于前轨迹l1的末尾，成为新的前轨迹l1的最后一点。重复该步骤，直到时间坐标t为丢失区域的中间帧。
[0122]
设j为位置预测参照帧数，取后轨迹l2的前j个点，对这j个点的时间坐标t和坐标x、y进行线性回归，获得坐标x、y相对于时间坐标t的线性方程：
[0123]
x(t)＝k
x
t b
x
[0124]
y(t)＝kyt by[0125]
其中k
x
为x相对于t的斜率，b
x
为x相对于t的截距，ky为y相对于t的斜率，by为y相对于t的截距。
[0126]
将线性方程x(t)、y(t)代入t＝f
21-1即得后轨迹l2向前一帧的线性预测位置的坐标x、y，其中f
21
为后轨迹l2的第一帧，这样获得的预测坐标记为(x
p
，y
p
)；
[0127]
在后轨迹l2向前一帧中，计算所有粒子的二维坐标(x
′
，y
′
)到预测坐标(x
p
，y
p
)的距离，距离(x
p
，y
p
)最近的(x
′
，y
′
)选为后轨迹l2向前一帧的粒子二维坐标，记为(xr，yr)，把(xr，yr)置于后轨迹l2的开头；
[0128]
若为三维轨迹则把(xr，yr，zf)置于后轨迹l2的开头，成为新的后轨迹l2的第一点，重复该步骤，直到时间坐标t为丢失区域的中间帧的后一帧。
[0129]
对l1的向前线性拟合帧和l2的向后线性拟合帧同样如此执行该步骤，如图7所示。
[0130]
第六步、若为三维轨迹则恢复拟合位置的三维信息。若样品图像含有三维信息，如全息图像、双目视觉、全内反射光强图等，则可在样品图像的(xr，yr)处，使用该三维信息对应的方式获取该粒子的z坐标，记为zr。则该粒子坐标为(xr，yr，zr)。
[0131]
以全息图像为例，通过高斯拟合获得粒子的三维坐标，设zg为高斯拟合范围，在坐标(xr，yr，zf)处，通过rayleigh-sommerfeld算法，对全息图像进行三维重建，获取z
f-zg到zf zg的三维光场，计算公式如下：
[0132][0133]
u(r，z)＝ft-1
(ft(is(r，0)
·
h(q，-z)))
[0134]
式中，h(r，-z)为瑞利索末菲传播子，r和z分别为拟合位置的横向坐标和纵向坐标，i为虚数单位，k为波数，r为光传播的距离，u(r，z)是三维重建所得的光场，is(r，0)为全息图像在r处的光强，ft-1
、ft分别为傅里叶逆变换和傅里叶变换，h(q，-z)为h(r，-z)的傅里叶变换，q为频域坐标。
[0135]
对光场中每一层不同的z，其中坐标(xr，yr)出对应的粒子所包含的所有像素的强度的和为iz，如图8所示。这样每个不同的z都有一个iz。对iz和z进行高斯拟合，拟合结果的
最大值记为zr，粒子的三维坐标即为(xr，yr，zr)，如图9所示。
[0136]
但第六步的实施方法并不限于全息图像，各种包含三维信息的图像均可应用，如双目视觉、全内反射光强图等。
[0137]
第七步，验证前轨迹l1和后轨迹l2是否属于同一个粒子，设纠错距离为wrong_d，计算前轨迹l1的最后一点到后轨迹l2的第一点的距离d
finish
，若大于wrong_d，则前轨迹l1和后轨迹l2不属于同一个粒子，否则前轨迹l1和后轨迹l2属于同一个粒子，如图10所示。
[0138]
续连前的轨迹碎片和续连完成的轨迹示例如图11所示。
[0139]
上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。