2. 立体匹配算法实现

任务目标

实现基于块匹配（Block Matching）的立体匹配算法，对双目图像对计算视差图，并定量评估算法在纹理丰富区域与弱纹理区域（如墙面、天空）的性能差异。

环境与依赖

Python 3.8+，依赖库：NumPy、OpenCV（含 contrib 模块）、Matplotlib（仅用于可视化验证）。

核心实现

以下代码实现固定窗口的块匹配视差计算，使用 SAD（差之和）作为相似度度量，并包含左-右一致性检查以剔除遮挡和误匹配点。视差搜索范围为 [0, max_disp)。

import numpy as np
import cv2

def block_matching_sad(left_gray, right_gray, max_disp=64, block_size=9):
    """
    基于SAD的块匹配立体匹配。
    返回: 视差图 (float32), 无效点置为 NaN; 以及有效像素掩码。
    """
    h, w = left_gray.shape
    pad = block_size // 2
    # 对右图进行零填充，保证窗口在边界处可计算
    right_padded = cv2.copyMakeBorder(right_gray, pad, pad, pad, pad,
                                       cv2.BORDER_CONSTANT, value=0)

    disp_map = np.full((h, w), np.nan, dtype=np.float32)
    valid_mask = np.zeros((h, w), dtype=bool)

    # 左图窗口块（避免重复提取）
    left_padded = cv2.copyMakeBorder(left_gray, pad, pad, pad, pad,
                                     cv2.BORDER_CONSTANT, value=0)

    for y in range(h):
        for x in range(w):
            # 提取左图窗口
            left_block = left_padded[y:y+block_size, x:x+block_size].astype(np.int16)

            best_cost = float('inf')
            best_d = 0
            # 在视差范围内搜索
            d_min = 0
            d_max = min(max_disp, x)  # 避免超出左图边界
            if d_max < d_min:
                continue

            for d in range(d_min, d_max + 1):
                # 右图对应窗口 (x - d)
                right_block = right_padded[y:y+block_size,
                                           x+pad-d : x+pad-d+block_size].astype(np.int16)
                # SAD 代价
                cost = np.sum(np.abs(left_block - right_block))
                if cost < best_cost:
                    best_cost = cost
                    best_d = d

            # 唯一性约束：次优代价应显著大于最优代价（此处简化为阈值判断）
            # 实际应用中可增加此检查，本例为保持核心简洁，仅做左右一致性
            disp_map[y, x] = float(best_d)
            valid_mask[y, x] = True

    # 左-右一致性检查（LRC）
    # 计算右图的视差图（通过交换左右图）
    disp_map_r = block_matching_sad_core(right_gray, left_gray, max_disp, block_size)

    # 一致性阈值 1 像素
    lr_diff = np.abs(disp_map - disp_map_r)
    consistent = (lr_diff <= 1.0) & np.isfinite(disp_map) & np.isfinite(disp_map_r)
    disp_map[~consistent] = np.nan
    valid_mask = consistent

    return disp_map, valid_mask

def block_matching_sad_core(left_gray, right_gray, max_disp, block_size):
    """内部函数：仅计算视差图，不做LRC，供LRC调用。"""
    h, w = left_gray.shape
    pad = block_size // 2
    right_padded = cv2.copyMakeBorder(right_gray, pad, pad, pad, pad,
                                       cv2.BORDER_CONSTANT, value=0)
    left_padded = cv2.copyMakeBorder(left_gray, pad, pad, pad, pad,
                                     cv2.BORDER_CONSTANT, value=0)
    disp_map = np.full((h, w), np.nan, dtype=np.float32)

    for y in range(h):
        for x in range(w):
            left_block = left_padded[y:y+block_size, x:x+block_size].astype(np.int16)
            best_cost = float('inf')
            best_d = 0
            d_max = min(max_disp, x)
            if d_max < 0:
                continue
            for d in range(0, d_max + 1):
                right_block = right_padded[y:y+block_size,
                                           x+pad-d : x+pad-d+block_size].astype(np.int16)
                cost = np.sum(np.abs(left_block - right_block))
                if cost < best_cost:
                    best_cost = cost
                    best_d = d
            disp_map[y, x] = float(best_d)
    return disp_map

def evaluate_texture_regions(disp_map, valid_mask, left_gray, texture_threshold=15.0):
    """
    基于左图梯度幅值区分纹理区域，评估视差图精度。
    返回: (纹理丰富区域有效比例, 弱纹理区域有效比例)
    """
    # 计算左图梯度幅值作为纹理强度
    grad_x = cv2.Sobel(left_gray, cv2.CV_64F, 1, 0, ksize=3)
    grad_y = cv2.Sobel(left_gray, cv2.CV_64F, 0, 1, ksize=3)
    grad_mag = np.sqrt(grad_x**2 + grad_y**2)

    texture_mask = grad_mag > texture_threshold
    weak_mask = ~texture_mask

    # 有效视差比例（非NaN且通过LRC）
    valid_disp = np.isfinite(disp_map) & valid_mask

    if np.sum(texture_mask) == 0 or np.sum(weak_mask) == 0:
        return 0.0, 0.0

    tex_valid_ratio = np.sum(valid_disp & texture_mask) / np.sum(texture_mask)
    weak_valid_ratio = np.sum(valid_disp & weak_mask) / np.sum(weak_mask)
    return tex_valid_ratio, weak_valid_ratio

# 最小调用示例（使用OpenCV示例图像）
if __name__ == "__main__":
    # 使用OpenCV自带的立体图像对（左/右）
    left = cv2.imread('left.png', cv2.IMREAD_GRAYSCALE)
    right = cv2.imread('right.png', cv2.IMREAD_GRAYSCALE)
    if left is None or right is None:
        # 若无图像，生成合成图像对用于验证（正弦纹理+平移）
        print("未找到图像，使用合成数据验证...")
        x = np.arange(320).reshape(1, -1)
        y = np.arange(240).reshape(-1, 1)
        pattern = (np.sin(x/10) + np.sin(y/15) * 2) * 30 + 128
        pattern = pattern.astype(np.uint8)
        # 模拟视差为8的平移
        left = pattern[:, 8:]
        right = pattern[:, :-8]
        # 填充至相同宽度
        left = cv2.copyMakeBorder(left, 0, 0, 0, 8, cv2.BORDER_REPLICATE)
        right = cv2.copyMakeBorder(right, 0, 0, 8, 0, cv2.BORDER_REPLICATE)

    disp, mask = block_matching_sad(left, right, max_disp=32, block_size=9)
    tex_ratio, weak_ratio = evaluate_texture_regions(disp, mask, left)

    print(f"纹理丰富区域有效视差比例: {tex_ratio:.2%}")
    print(f"弱纹理区域有效视差比例: {weak_ratio:.2%}")

    # 可视化（可选）
    import matplotlib.pyplot as plt
    plt.figure(figsize=(12, 4))
    plt.subplot(131), plt.imshow(left, cmap='gray'), plt.title('Left')
    plt.subplot(132), plt.imshow(disp, cmap='jet'), plt.title('Disparity')
    plt.subplot(133), plt.imshow(mask, cmap='gray'), plt.title('Valid Mask')
    plt.tight_layout()
    plt.savefig('bm_result.png', dpi=100)
    print("结果图已保存为 bm_result.png")

运行与验证

运行方式：将代码保存为 block_matching.py，在项目目录放置双目图像 left.png 和 right.png（校正后、行对齐），执行 python block_matching.py。若无测试图像，程序自动生成具有已知视差（8像素）的合成正弦纹理图像对进行验证。

验证方法：

合成数据验证：合成图像对具有真实视差8像素，程序输出视差图中应大部分区域视差值接近8，且有效比例接近100%。
真实图像评估：观察输出的 bm_result.png，视差图中物体边缘应清晰，平坦区域（弱纹理）会出现较多黑色空洞（无效点）。控制台输出的两个比例应满足：纹理丰富区域有效比例显著高于弱纹理区域（通常高30个百分点以上），这验证了块匹配在纹理丰富处性能更优。
判定通过条件：合成数据测试中，有效像素中视差误差≤1像素的比例超过90%；真实图像评估中，纹理丰富区域有效比例 > 60% 且弱纹理区域有效比例 < 40%。

常见问题

运行缓慢：纯Python双循环在较大图像上耗时较长。可先用 cv2.resize 缩小图像测试，或改用C++实现。
视差图出现条纹：通常因窗口过小或图像噪声导致，可增大 block_size（如15）或先做高斯滤波。
左右一致性检查后有效点过少：可放宽阈值至2像素，或检查输入图像是否已校正（极线对齐）。

📚 课程目录