细菌细胞追踪：帧间关联算法详解与代码实现

2025-03-03 04:06:59

追踪细菌细胞轮廓：解决帧间关联问题

这个问题挺常见的：你有连续两帧细菌细胞的图像，想计算每个细胞的瞬时速度。你知道每帧图像里每个细胞的质心 (COM) 坐标，但不知道如何把两帧图像中对应的细胞关联起来，也就是，第一帧的哪个细胞是第二帧的哪个细胞。

我理解你尝试过用轮廓特征（如长短轴）给每个细胞一个独一无二的 ID。这思路不错，但前提是图像分割非常完美，且每个细菌细胞形态差异足够大。实际情况往往不是这样，对吧？所以我们需要更好的办法。

为什么直接匹配特征不可靠？

分割不完美： 你的代码找到的轮廓可能不准，一个细菌可能被分成多个轮廓，或者多个细菌被识别成一个。
细胞相似： 细菌长得都差不多，大小、形状都可能很接近，仅靠这些特征很难区分。
细胞运动： 细菌会移动、旋转，甚至形状会略微改变，导致特征变化。

可行的解决方案

下面几个方法，你可以试试，它们各有优缺点，可以根据实际情况组合使用：

1. 最近邻匹配 (Nearest Neighbor Matching)

原理： 这是最简单粗暴的办法。对第一帧的每个细胞，在第二帧里找离它最近的细胞，认为它们是同一个。
代码示例：

import numpy as np
from scipy.spatial import distance

def nearest_neighbor_matching(centroids1, centroids2):
    """
    使用最近邻算法匹配两组质心。

    参数：
        centroids1: 第一帧的质心坐标，形状为 (N1, 2) 的 NumPy 数组。
        centroids2: 第二帧的质心坐标，形状为 (N2, 2) 的 NumPy 数组。

    返回值：
        matches: 一个列表，长度为 N1。
                 matches[i] 是第二帧中与第一帧第 i 个质心匹配的质心的索引。
                 如果找不到匹配项，则为 -1。
    """
    matches = []
    for i, centroid1 in enumerate(centroids1):
        distances = distance.cdist([centroid1], centroids2) # 计算距离
        nearest_index = np.argmin(distances)
        matches.append(nearest_index)
    return matches

# 假设你已经有了两帧图像的质心坐标：
# centroids1 = [(x1_1, y1_1), (x1_2, y1_2), ...]  # 第一帧的质心
# centroids2 = [(x2_1, y2_1), (x2_2, y2_2), ...]  # 第二帧的质心

# 使用示例：
centroids1 = np.array([[10, 20], [30, 40], [50, 60]]) #示例数据
centroids2 = np.array([[12, 22], [35, 38], [53, 61]])
matches = nearest_neighbor_matching(centroids1, centroids2)
print(matches) # 示例: 可能输出 [0, 1, 2], 表示第一帧的每个细胞，和第二帧的对应序号细胞相匹配
#得到匹配序号后, 便可以着手计算速度等后续操作

操作步骤：
1. 提取两帧图像中所有细胞的质心坐标。
2. 用 nearest_neighbor_matching 函数进行匹配。
3. matches列表给出的就是匹配关系，matches[i]对应第二帧中的细胞编号。
优点： 简单、快速。
缺点： 如果细胞运动距离大于细胞间距，或者细胞密度太大，很容易出错。
安全建议 :
- 在计算速度之前, 可对匹配结果做过滤。如设定距离阈值, 当两帧质心距离过大时, 视为匹配无效, 不参与后续计算。

2. 基于匈牙利算法的匹配 (Hungarian Algorithm)

原理： 匈牙利算法（也叫 Kuhn-Munkres 算法）能解决“分配问题”。想象一下，你有 N 个任务，N 个工人，每个工人做每个任务的成本不一样。匈牙利算法能找到总成本最低的任务分配方案。在这里，第一帧的细胞是“任务”，第二帧的细胞是“工人”，“成本”就是细胞之间的距离。
代码示例：

import numpy as np
from scipy.optimize import linear_sum_assignment
from scipy.spatial import distance

def hungarian_matching(centroids1, centroids2):
    """
    使用匈牙利算法匹配两组质心。

    参数：
        centroids1: 第一帧的质心坐标，形状为 (N1, 2) 的 NumPy 数组。
        centroids2: 第二帧的质心坐标，形状为 (N2, 2) 的 NumPy 数组。

    返回值：
        matches: 一个列表，长度为 min(N1, N2)。
                 每个元素是一个元组 (i, j)，表示第一帧的第 i 个质心与第二帧的第 j 个质心匹配。
    """

    # 计算成本矩阵（距离矩阵）
    cost_matrix = distance.cdist(centroids1, centroids2)

    # 使用匈牙利算法
    row_ind, col_ind = linear_sum_assignment(cost_matrix)

    # 构建匹配结果
    matches = list(zip(row_ind, col_ind))
    return matches

# 假设你已经有了两帧图像的质心坐标：
# centroids1 = [(x1_1, y1_1), (x1_2, y1_2), ...]  # 第一帧的质心
# centroids2 = [(x2_1, y2_1), (x2_2, y2_2), ...]  # 第二帧的质心
# 使用示例：
centroids1 = np.array([[10, 20], [30, 40], [50, 60],[99,99]]) # 示例数据
centroids2 = np.array([[12, 22], [35, 38], [53, 61]])

matches = hungarian_matching(centroids1, centroids2)
print(matches)  # 示例: 可能输出 [(0, 0), (1, 1), (2, 2)],
                #表示第一帧的第 0,1,2 号细胞和第二帧的 0,1,2 号细胞相匹配。多出的[99,99]因为无匹配对象而被舍弃。

操作步骤：
1. 提取两帧图像中所有细胞的质心坐标。
2. 用 hungarian_matching 函数进行匹配。
3. matches 列表给出的就是匹配关系。
优点： 比最近邻匹配更准确，尤其在细胞密度较大时。
缺点： 如果两帧中细胞数量不同，有一部分细胞会无法匹配。计算量比最近邻大。

3. 卡尔曼滤波 (Kalman Filter)

原理： 卡尔曼滤波是一种强大的追踪算法。它不仅考虑当前帧的信息，还会利用之前的运动状态来预测当前位置，所以对噪声和遮挡有一定鲁棒性。每个细胞都有一个卡尔曼滤波器，跟踪它的位置、速度等。
代码示例： (因为涉及到状态方程,代码较为复杂，这里只给出简化的、使用 filterpy 库的示例)：

import numpy as np
from filterpy.kalman import KalmanFilter
from filterpy.common import Q_discrete_white_noise
from scipy.spatial import distance

def initialize_kalman_filter(initial_position):
    """
    初始化一个卡尔曼滤波器。

    参数：
        initial_position: 细胞的初始位置 (x, y)。

    返回值：
        kf: 一个 KalmanFilter 对象。
    """

    # 状态向量：[x, y, vx, vy] (位置和速度)
    kf = KalmanFilter(dim_x=4, dim_z=2)

    # 初始状态
    kf.x = np.array([initial_position[0], initial_position[1], 0, 0])

    # 状态转移矩阵 (假设匀速运动模型)
    dt = 1  # 时间间隔，可以根据你的视频帧率调整
    kf.F = np.array([[1, 0, dt, 0],
                     [0, 1, 0, dt],
                     [0, 0, 1, 0],
                     [0, 0, 0, 1]])

    # 测量矩阵 (我们直接测量位置)
    kf.H = np.array([[1, 0, 0, 0],
                     [0, 1, 0, 0]])

    # 测量噪声协方差 (可以根据实际情况调整)
    kf.R = np.array([[5, 0],
                     [0, 5]])

    # 过程噪声协方差 (可以根据实际情况调整)
    kf.Q = Q_discrete_white_noise(dim=4, dt=dt, var=0.1)

    # 初始状态协方差 (可以根据实际情况调整)
    kf.P = np.eye(4) * 10

    return kf

def kalman_tracking(centroids_list):
    """
   使用卡尔曼滤波进行多目标追踪。

    参数：
        centroids_list: 每一帧的质心列表。例如：
                        [[centroid1_frame1, centroid2_frame1, ...],
                         [centroid1_frame2, centroid2_frame2, ...],
                         ...]

    返回值：
        tracks: 一个字典，键是细胞的 ID，值是每个细胞的轨迹（一系列位置）。
    """
    trackers = {}
    track_id = 0
    tracks = {}
    #遍历每一帧
    for frame_num, centroids in enumerate(centroids_list):
        if frame_num == 0: #第0帧，初始化追踪器
            #对第一个frame做一些处理
            for centroid in centroids:
                trackers[track_id] = initialize_kalman_filter(centroid) #为这个点创建一个追踪器。
                tracks[track_id] = [centroid] #把第一个点加进track
                track_id+=1

        else: #在第一帧之后
            # 先进行预测
            for kf in trackers.values():
                kf.predict()

            #用这一帧的所有观测点，更新所有的追踪器。
            #但在这之前我们需要根据预测结果和观测结果计算匹配关系.

            # 成本矩阵: 预测位置 与 当前帧观测的质心 的距离
            predicted_positions = np.array([kf.x[:2] for kf in trackers.values()])
            if(len(predicted_positions) >0 and len(centroids)>0):
                cost_matrix = distance.cdist(predicted_positions, centroids)

                # 使用匈牙利算法进行匹配 (与上面的方法 2 相同)
                row_ind, col_ind = linear_sum_assignment(cost_matrix)
                matches = list(zip(row_ind, col_ind))
            else:
                matches=[] #没有匹配目标

            # 记录未匹配的
            unmatched_trackers = set(trackers.keys()) - {row_ind[i] for i in range(len(row_ind))}
            unmatched_centroids = set(range(len(centroids))) - {col_ind[i] for i in range(len(col_ind))}

             # 更新已匹配的追踪器
            for tracker_index, centroid_index in matches:
                tracker_id = list(trackers.keys())[tracker_index]
                trackers[tracker_id].update(centroids[centroid_index])
                tracks[tracker_id].append(centroids[centroid_index])

              # 对未匹配到的 新观测点：添加新追踪器。
            for i in unmatched_centroids:
                trackers[track_id] = initialize_kalman_filter(centroids[i])
                tracks[track_id] = [centroids[i]]
                track_id += 1

             #对于未匹配到的追踪器, 可做额外处理
            for tracker_id in unmatched_trackers:
                #pass #简单忽略.
                tracks[tracker_id].append(None)  # 可用None,表示此帧中丢了。


    return tracks

# 示例：(注意, 数据需提前准备好。)
centroids_frame1 = np.array([[10, 20], [30, 40], [50, 60]])
centroids_frame2 = np.array([[12, 22], [31, 41], [53, 59],[99,99]]) #新增的细菌。
centroids_frame3 = np.array([[13, 24], [33, 40], [55, 62]])#99消失
#centroids_list=[centroids_frame1] #仅用于测试匈牙利算法时使用。
centroids_list = [centroids_frame1,centroids_frame2,centroids_frame3] # 注意这是所有帧
tracks = kalman_tracking(centroids_list)

# 现在, 'tracks'里应该包含了每个细胞的轨迹。
print (tracks)

for track_id, track in tracks.items():
    print(f"细胞 {track_id}: {track}")
#可以输出类似:
# 细胞 0: [array([10, 20]), array([12, 22]), array([13, 24])]
#这样的结果。

操作步骤：
1. 对第一帧的每个细胞，初始化一个卡尔曼滤波器。
2. 对后续每一帧：
  - 用卡尔曼滤波器预测每个细胞的位置。
  - 用匈牙利算法将预测位置与当前帧的观测位置（质心）进行匹配。
  - 用观测位置更新卡尔曼滤波器的状态。
  - 对未匹配的观测位置,认为可能是新出现的细胞，添加新的卡尔曼滤波器。
3. 不断迭代，每个细胞都会被持续追踪
进阶技巧：
- 非线性运动模型： 如果细胞运动不规律，可以考虑更复杂的运动模型，比如“恒定加速度模型”或“交互多模型”(IMM)。
- 数据关联： 除了匈牙利算法，还可以尝试更高级的数据关联算法，比如“联合概率数据关联”(JPDA) 或“多假设跟踪”(MHT)。
- 外观信息： 除了位置，还可以把细胞的形状、大小等特征融入卡尔曼滤波，提高追踪准确性。可以修改测量向量和测量矩阵。
安全建议
- 调整kf.R和kf.Q, 这是调参重点, 相信可以做出比最近邻算法更好的结果.
- 注意数据溢出和除零错误.

处理图像预处理的代码建议:

你的预处理代码有些地方值得改进. 可以让结果更好：
1. gaussian_filter_multiscale_retinex 这个函数：
* 返回值可以不用归一化到 0-255，直接用 float 类型，保留更多信息。

2. `process_image` 这个函数：
   * 可以考虑使用形态学操作（开运算、闭运算）来去除噪点、连接断开的轮廓。
   * `cv.fitEllipse` 可能有时会失效 (返回奇怪的值), 导致`calculate_ellipse_area`出错。建议增加异常值检测与处理
   * 可以考虑用 `cv.contourArea` 直接计算轮廓面积, 更可靠，快速。