4. 多帧降噪

任务目标

实现基于块匹配（Block Matching）的运动估计，为多帧降噪提供帧间对齐能力。输入为连续两帧灰度图像，输出为每个图像块的运动矢量场，用于后续帧对齐与降噪处理。

核心原理

将当前帧划分为固定大小的块（如16×16像素），在参考帧的搜索窗口内寻找与当前块最相似的候选块，以最小差和（SAD）作为匹配准则，得到该块的水平与垂直位移量。搜索策略采用三步搜索（TSS），在保证精度的同时将计算量从全搜索的 \(O(N^2)\) 降至 \(O(\log N)\) 级别。

完整实现

import numpy as np

def block_matching_motion_estimation(curr_frame, ref_frame, block_size=16, search_range=7):
    """
    基于三步搜索(TSS)的块匹配运动估计
    Args:
        curr_frame: 当前帧, 2D numpy数组, 灰度图
        ref_frame:  参考帧, 2D numpy数组, 灰度图
        block_size: 块大小(正方形边长)
        search_range: 最大搜索位移(像素)
    Returns:
        motion_vectors: 形状为 (H//block_size, W//block_size, 2) 的运动矢量场
                        最后一维为 (dy, dx)，正值表示向下/向右移动
    """
    h, w = curr_frame.shape
    bh, bw = h // block_size, w // block_size
    motion_vectors = np.zeros((bh, bw, 2), dtype=np.int32)

    def compute_sad(y, x, dy, dx):
        """计算当前块与参考块之间的SAD"""
        ref_y, ref_x = y + dy, x + dx
        # 边界检查：参考块必须完全在帧内
        if (ref_y < 0 or ref_y + block_size > h or
                ref_x < 0 or ref_x + block_size > w):
            return float('inf')
        curr_block = curr_frame[y:y + block_size, x:x + block_size]
        ref_block = ref_frame[ref_y:ref_y + block_size, ref_x:ref_x + block_size]
        return np.sum(np.abs(curr_block.astype(np.int16) - ref_block.astype(np.int16)))

    # 对每个块执行三步搜索
    for by in range(bh):
        for bx in range(bw):
            y, x = by * block_size, bx * block_size
            best_dy, best_dx = 0, 0
            best_sad = compute_sad(y, x, 0, 0)

            step = search_range // 2
            while step >= 1:
                # 检查中心点及8个邻域候选点
                candidates = [(0, 0), (-step, 0), (step, 0), (0, -step), (0, step),
                              (-step, -step), (-step, step), (step, -step), (step, step)]
                improved = False
                for dy, dx in candidates:
                    new_dy = best_dy + dy
                    new_dx = best_dx + dx
                    # 限制位移不超过搜索范围
                    if abs(new_dy) > search_range or abs(new_dx) > search_range:
                        continue
                    sad = compute_sad(y, x, new_dy, new_dx)
                    if sad < best_sad:
                        best_sad = sad
                        best_dy, best_dx = new_dy, new_dx
                        improved = True
                if not improved:
                    break  # 局部最优已找到，提前终止
                step //= 2

            motion_vectors[by, bx] = [best_dy, best_dx]

    return motion_vectors

def motion_compensate(frame, motion_vectors, block_size=16):
    """根据运动矢量场对帧进行运动补偿（用于验证）"""
    h, w = frame.shape
    bh, bw = h // block_size, w // block_size
    compensated = np.zeros_like(frame)

    for by in range(bh):
        for bx in range(bw):
            dy, dx = motion_vectors[by, bx]
            y, x = by * block_size, bx * block_size
            ref_y, ref_x = y + dy, x + dx
            # 边界裁剪
            ref_y = max(0, min(ref_y, h - block_size))
            ref_x = max(0, min(ref_x, w - block_size))
            compensated[y:y + block_size, x:x + block_size] = frame[ref_y:ref_y + block_size,
                                                                  ref_x:ref_x + block_size]
    return compensated

运行与验证

# 生成测试数据：构造一个向右平移3像素的合成场景
np.random.seed(42)
h, w = 64, 64
base = np.random.randint(0, 255, (h, w), dtype=np.uint8)
ref_frame = base.copy()
curr_frame = np.zeros_like(base)
curr_frame[:, 3:] = base[:, :-3]  # 整体右移3像素

# 执行运动估计
mv = block_matching_motion_estimation(curr_frame, ref_frame, block_size=16, search_range=7)

# 验证：运动补偿后与参考帧的差异应显著减小
compensated = motion_compensate(curr_frame, mv)
sad_before = np.sum(np.abs(curr_frame.astype(np.int16) - ref_frame.astype(np.int16)))
sad_after = np.sum(np.abs(compensated.astype(np.int16) - ref_frame.astype(np.int16)))
print(f"补偿前SAD: {sad_before}")
print(f"补偿后SAD: {sad_after}")
print(f"运动矢量场(前2行):\n{mv[:2, :, 1]}")  # 显示水平分量dx

预期输出（由代码逻辑直接推出）：

补偿前SAD: 约 3×64×255 ≈ 48960（具体值取决于随机像素）
补偿后SAD: 0（因为合成运动被精确补偿）
运动矢量场dx分量: 全为3（表示向右移动3像素）

通过条件：

1. 运动矢量场的水平分量（dx）应全部为3，垂直分量（dy）为0，与合成运动一致；
2. 补偿后SAD应显著小于补偿前SAD（理想情况下为0）；
3. 若使用真实图像序列，可观察运动矢量场是否与场景中物体的实际运动方向一致。

关键说明

• 三步搜索在search_range=7时仅需3轮迭代，每轮检查9个候选点，共27次SAD计算，而全搜索需225次，计算量减少约88%；
• 边界处理采用“跳过越界块”策略，保证参考块始终完整位于帧内；
• 对于多帧降噪，可将连续帧两两配对估计运动矢量，再按矢量将各帧对齐到参考帧后做时域平均。