COLMAP PatchMatch Stereo 算法详解

COLMAP 的 PatchMatch Stereo 是 MVS 稠密重建的核心算法，通过随机初始化深度和法线、NCC 代价评估、空间传播与随机扰动迭代优化，估计每个像素的最优深度。本文记录算法原理与实现细节。

算法核心思想

用一句话概括：每个像素随机猜一个深度，靠 NCC 判断猜得对不对，好的猜测向邻居扩散，反复多轮直到收敛。

展开来说有三个机制：

1. 随机猜测（探索）

每个像素初始时在 [depth_min, depth_max] 范围内随机取一个深度值——纯粹蒙的。后续每轮迭代也在当前最优值附近随机扰动，探索新的可能性。

2. NCC 打分（判断对错）

判断一个深度猜测是否正确的依据：假设该像素深度为 d，算出 3D 点，投影到源图，比较参考图 patch 和源图 patch：

深度正确 → 投影到源图是同一块表面 → patch 相似 → NCC 高 → 猜测"好"
深度错误 → 投影到源图是不相关的区域 → patch 不像 → NCC 低 → 猜测"差"

3. 空间传播（好的答案扩散）

如果邻居像素的深度比我好（NCC 更高），我就直接拿来用。通过多轮迭代 + 四方向扫描，正确深度从少数碰巧猜对的像素逐步扩散到全图：

初始:  全图随机猜 → 几乎全错
  ↓
传播+打分循环:  好的猜测被 NCC 识别，通过传播扩散给邻居
  ↓
收敛:  大部分像素找到了 NCC 最高的深度

这正是 PatchMatch 名称的由来——通过 Patch 匹配（NCC）来驱动 Match 传播（空间扩散）。

算法定位

COLMAP MVS 流水线:

  image_undistorter          ← 去畸变，统一为 pinhole 模型
        │
  PatchMatch Stereo          ← 本文档重点
    ├── Phase 1: Photometric ← 仅光度一致性 (NCC)
    └── Phase 2: Geometric   ← 加入几何一致性约束
        │
  stereo_fusion              ← 多视角深度图融合
        │
  Poisson/Delaunay 网格化    ← (可选) 生成 mesh

输入与输出

输入	说明
参考图像 (ref)	待估计深度的图像
源图像集合 (src)	与参考图像有视角重叠的邻近图像（默认选取 20 张）
相机内参 K	每幅图像的焦距、主点
相对位姿 R, T	参考图像到每张源图像的旋转和平移
depth_min / depth_max	深度搜索范围

输出	说明
depth_map	每像素的深度值 (沿光轴距离)
normal_map	每像素的法线向量 (3D 单位向量)
consistency_graph	(可选) 多视角一致性图

可配置参数

源码位置: patch_match_options.h

参数	默认值	说明
max_image_size	-1	图像最大边长，超出则缩放
window_radius	5	NCC 窗口半径，实际窗口 = 2×5+1 = 11×11
window_step	1	NCC 窗口内像素步长
sigma_spatial	-1	双边空间权重 σ（默认 = window_radius）
sigma_color	0.2	双边颜色权重 σ
ncc_sigma	0.6	NCC 似然函数的扩散参数
num_samples	15	每像素每轮 Monte Carlo 源图采样次数
num_iterations	5	外层迭代次数（每迭代 = 4 个 sweep 方向）
min_triangulation_angle	1.0°	源图选择的最小三角化角度阈值
incident_angle_sigma	0.9	入射角似然函数的扩散参数
depth_min / depth_max	-1.0	深度搜索范围（从 sparse SfM 估计）
filter_min_ncc	0.1	过滤：最小 NCC 阈值
filter_min_triangulation_angle	3.0°	过滤：最小三角化角度
filter_min_num_consistent	2	过滤：最少一致源图数
geom_consistency	true	是否启用几何一致性
geom_consistency_regularizer	0.3	几何一致性项权重
geom_consistency_max_cost	3.0 px	最大前向-后向重投影误差
filter	true	是否启用输出过滤
filter_geom_consistency_max_cost	1.0 px	过滤时的几何一致性重投影误差阈值
cache_size	32.0 GB	工作空间缓存大小
gpu_index	“-1”	GPU 索引，“-1” 使用所有可用 GPU
write_consistency_graph	false	是否输出一致性图
allow_missing_files	false	是否容忍缺失文件

窗口半径上限 32（硬编码 kMaxPatchMatchWindowRadius），与 CUDA 线程块大小 THREADS_PER_BLOCK = 32 匹配。

数据结构

GPU 全局状态

源码位置: patch_match_cuda.h

变量	类型	尺寸	用途
ref_image_->image	GpuMat<uint8_t>	W×H	参考图灰度
ref_image_->sum_image	GpuMat<float>	W×H	预计算的 bilateral 加权和
ref_image_->squared_sum_image	GpuMat<float>	W×H	预计算的 bilateral 加权平方和
depth_map_	GpuMat<float>	W×H	每像素深度假设
normal_map_	GpuMat<float>	W×H×3	每像素法线假设 (nx, ny, nz)
cost_map_	GpuMat<float>	W×H×N_src	每源图的光度代价 (1-NCC)
sel_prob_map_	GpuMat<float>	W×H×N_src	当前源图选择概率
prev_sel_prob_map_	GpuMat<float>	W×H×N_src	上一 sweep 方向的选择概率
consistency_mask_	GpuMat<uint8_t>	W×H×N_src	一致性过滤掩码（仅最后一轮 sweep 分配）
src_depth_maps_texture_	Layered Texture (float)	源图深度图纹理（几何一致性模式启用时使用）
rand_state_map_	GpuMatPRNG	W×H	每像素 CUDA 随机数状态 (curandState)
global_workspace_	GpuMat<float>	max(W,H)×N_src×2	前向消息 + 采样 CDF

纹理缓存

源码位置: cuda_texture.h

变量	类型	说明
ref_image_texture_	Layered Texture (uint8)	参考图像，支持硬件双线性插值
src_images_texture_	Layered Texture (uint8)	所有源图像 (max_W × max_H × N_src)
poses_texture_[4]	Layered Texture (float)	4 种旋转的相机位姿，每种存储 N_src 个位姿

每个位姿存储 43 个 float：

[0..3]    : 源图内参 K_src (fx, cx, fy, cy)
[4..12]   : 相对旋转 R (3×3)
[13..15]  : 相对平移 T (3)
[16..18]  : 投影中心 C (3)
[19..30]  : 投影矩阵 P = K_src [R|T] (3×4)
[31..42]  : 逆投影矩阵 P⁻¹ (3×4)

初始化阶段

源图像选择

源码位置: patch_match.cc:291-368

从 patch-match.cfg 文件中读取参考图像和候选源图的配置。支持三种模式：

1. __all__ — 所有其他图像作为源图
2. __auto__, N — 自动选择最优 N 张源图（默认 N=20）
   • 按共享稀疏点数量排序
   • 过滤条件：第 75 百分位三角化角 > min_triangulation_angle
3. 显式列表 — 手动指定源图名称

双边预过滤

源码位置: gpu_mat_ref_image.cu:39-81

在迭代开始前，对参考图像的每个像素预计算 bilateral 加权的局部统计量：

for dr in range(-radius, radius+1, step):
    for dc in range(-radius, radius+1, step):
        ref_color = ref_image[row+dr, col+dc]

        # 双边权重: 空间距离 + 颜色差异
        w = exp(-(dr²+dc²) / (2σ_s²) - (center_color - ref_color)² / (2σ_c²))

        color_sum        += w * ref_color
        color_squared_sum += w * ref_color²
        weight_sum       += w

sum_image[row, col]        = color_sum / weight_sum
squared_sum_image[row, col] = color_squared_sum / weight_sum

目的：将参考图像侧的 NCC 统计量预计算好，迭代时只需计算源图像侧，大幅减少重复计算。

为什么参考图侧可以预计算：NCC 公式展开后包含 5 项累加，其中参考图侧的 2 项在整个迭代过程中始终不变——不管 depth/normal 怎么变，参考图上每个像素的颜色值是固定的，变化的只有 warp 到源图后读出的颜色：

NCC 的 5 项累加:

  参考图侧（不随假设变化，可预计算）:
    ref_sum     = Σ w(i) × ref_color(i)       ← 存入 sum_image
    ref_sq_sum  = Σ w(i) × ref_color(i)²      ← 存入 squared_sum_image

  源图侧（每个假设都要重新算）:
    src_sum      = Σ w(i) × src_color(i)      ← warp 改变 → 颜色改变
    src_sq_sum   = Σ w(i) × src_color(i)²     ← 同上
    cross_sum    = Σ w(i) × ref(i) × src(i)   ← 同上

NCC 窗口有 (2×5+1)² = 121 个像素。每个假设评估时，预计算省掉了 2 项 × 121 像素的累加，只需从全局内存读取 2 个 float（sum_image[row,col] 和 squared_sum_image[row,col]）。整个算法每个像素评估 5假设 × 15采样 × 20次迭代 = 1500 次假设，累计节省约 40% 的 NCC 计算量。

深度与法线初始化

深度 (Photometric 模式)：在 [depth_min, depth_max] 范围内均匀随机采样。

法线：使用 Marsaglia 方法在面向相机的半球面上均匀随机采样：

1. 采样 u₁, u₂ ∈ [-1, 1]，直到 u₁² + u₂² < 1
2. 法线 n = (2u₁√(1-u₁²-u₂²), 2u₂√(1-u₁²-u₂²), 1-2(u₁²+u₂²))
3. 确保法线朝向相机：dot(n, view_ray) < 0

选择概率：所有源图初始化为均匀概率 0.5。

核心算法：Sweep Kernel

源码位置: patch_match_cuda.cu:899-1254

整体调度

ComputeInitialCost()              ← 评估初始随机 depth/normal 的 NCC

for iter = 0 to num_iterations-1:   ← 默认 5 次迭代
    for sweep = 0 to 3:              ← 4 个方向
        perturbation = 1.0 / 2^(iter + sweep/4)
        prev_weight = (iter*4 + sweep) / (num_iterations*4)
        SweepFromTopToBottom(...)
        Rotate()                      ← 所有数据结构旋转 90°

每次 Rotate() 将整个 GPU 状态旋转 90°（depth_map、normal_map、cost_map、参考图像等），同时切换到预计算的旋转后位姿。这样 4 个 sweep 分别对应：上→下、右→左、下→上、左→右。

深度估计的基本原理

整个算法的核心问题是：怎么从多张照片判断一个像素的深度？

人眼通过左右眼的视差感知深度——两眼看到的画面有细微差异，差异越大说明物体越近。MVS 做同样的事，用"参考图 + N 张源图"代替"左眼 + 右眼"。判断流程为：

第 1 步: 对参考图像素 (row, col)，假设一个深度 d，算出 3D 点

    像素 (500, 300), 深度假设 d = 10m → 3D 点 P = (2.3, 1.8, 10.0)

第 2 步: 把 3D 点 P 投影到源图（通过相机位姿）

    源图像素 (480, 290)  ← 位置不同，因为视角不同

第 3 步: 比较参考图 patch 和源图 patch

    参考图 (500,300) 附近 11×11 区域 vs 源图 (480,290) 附近 11×11 区域

    d 正确: 两个 patch 看到同一块表面 → NCC ≈ 1（好）
    d 错误: 投影到了源图中不相关的位置 → NCC ≈ 0（差）

猜错了的例子：

真实情况: 像素 (500,300) 对应一堵白墙，深度 10m

假设 d = 10m (正确):
  3D 点在白墙上 → 投影到源图也是白墙 → patch 相似 → NCC = 0.95 ✓

假设 d = 5m (太近):
  3D 点在白墙前面 → 投影到源图时落到墙前的一棵树上
  → 参考图 patch 是白墙，源图 patch 是树叶 → NCC = 0.1 ✗

假设 d = 20m (太远):
  3D 点穿过了墙 → 投影到源图时落到后面的建筑上
  → 参考图 patch 是白墙，源图 patch 是建筑 → NCC = 0.2 ✗

判断"好还是不好"的标准就是 NCC——深度正确时，参考图和源图的对应 patch 看起来一样（高 NCC）；深度错误时，投影跑偏了，patch 不匹配（低 NCC）。整个算法就是反复猜不同的 depth + normal，用 NCC 给每个猜测打分，留下得分最高的。

Sweep 是什么

一次 sweep 就是从一端扫到另一端，逐像素处理。以从上到下的 sweep 为例：

第  0 行: 处理所有像素 (更新深度)
第  1 行: 处理所有像素 (可以看第 0 行的结果)
第  2 行: 处理所有像素 (可以看第 1 行的结果)
...
最后一行: 处理所有像素

扫一遍的过程中，每个像素都做了"看看上面邻居的深度好不好，好就拿来用"的操作。

为什么要多轮迭代

初始状态每个像素的深度都是随机猜的，全不对。靠传播让正确值从第 0 行逐行接力传到第 1536 行需要 1536 次 sweep，根本等不起。

PatchMatch 的关键是不只依赖传播，还同时做随机扰动：

每轮 sweep 做两件事:

1. 传播: 抄邻居的好结果 → 正确信息向前推进一步
2. 随机扰动: 每个像素自己也在探索新深度 → 可能直接猜中正确答案

两个机制叠加:
  传播让正确值稳步扩散（确定性推进）
  随机扰动让远处像素有机会自己找到正确答案（不需要等传播传过来）

类比一个教室里老师把答案告诉第一排，每个学生可以抄前一排的答案（传播），也可以自己猜（随机扰动），比较后留下更好的。几轮之后全班都有正确答案——不需要等第一排的答案逐排传到最后一排。

所以 5 轮迭代 × 4 方向 = 20 次 sweep 就够了：

初始:  全部随机 → 全错
第 1 轮:  少数像素碰巧猜对 + 传播到邻居 → 只覆盖少数区域
第 2 轮:  更多像素自己猜对 + 传播继续扩散 → 覆盖更多区域
...
第 5 轮:  大部分像素已收敛，只剩微调

为什么要 4 个方向扫描

只从上往下扫，信息只能向下传播，下方的正确值永远传不到上方。所以扫 4 个方向，让信息从四面八方传播：

sweep 0: 上→下   （上方好结果往下传）
sweep 1: 右→左   （右方好结果往左传）
sweep 2: 下→上   （下方好结果往上传）
sweep 3: 左→右   （左方好结果往右传）

4 个方向扫完 = 信息在上下左右都传播了一轮。5 次迭代 × 4 方向 = 20 轮扫描。

为什么用"旋转 90°"实现 4 方向

最直觉的做法是写 4 个不同的 kernel：SweepDown()、SweepRight()、SweepUp()、SweepLeft()。

COLMAP 只写了 1 个 kernel（SweepFromTopToBottom），永远只从上往下扫。每扫完一次把整张图旋转 90°——kernel 还是同一个，但效果等价于换了方向扫原图：

原始图:           旋转 90° 后:      再旋转 90°:      再旋转 90°:
┌──────┐         ┌──────┐        ┌──────┐        ┌──────┐
│↓ ↓ ↓ │         │↓ ↓ ↓ │        │↓ ↓ ↓ │        │↓ ↓ ↓ │
│↓ ↓ ↓ │  →      │↓ ↓ ↓ │  →     │↓ ↓ ↓ │  →     │↓ ↓ ↓ │
│↓ ↓ ↓ │         │↓ ↓ ↓ │        │↓ ↓ ↓ │        │↓ ↓ ↓ │
└──────┘         └──────┘        └──────┘        └──────┘
上→下            上→下            上→下            上→下
(实际:上→下)    (实际:右→左)    (实际:下→上)    (实际:左→右)

好处：只需写、调试、维护 1 个 kernel 而非 4 个。

perturbation 和 prev_weight

• perturbation：每轮随机扰动的幅度。第 1 轮大胆猜（搜索范围 ±100%），第 5 轮只微调（±6%）。从粗到精。
• prev_weight：选择源图时多大程度信任上一轮的结果。早期 trust 当前计算，后期 trust 历史积累。防止结果抖动。

线程组织

• 每 CUDA 线程处理一列，从上到下逐行扫描
• 32 线程/块 (THREADS_PER_BLOCK = 32)
• 同一块内的线程共享参考图像 tile（共享内存）
• 第一行加载完整 tile，后续行增量更新（移出旧行，加载新行）

Sweep 内逐像素处理流程

每个像素 (row, col) 的处理分为以下步骤：

步骤 ①：后向消息预计算

源码位置: patch_match_cuda.cu:942-955

在 sweep 开始前，为每张源图从下到上计算后向消息 β：

对每张源图 image_idx:
    beta = 0.5   (均匀初始)
    for row = H-1 down to 0:
        cost = cost_map[row, col, image_idx]
        beta = ComputeBackwardMessage(cost, beta)
        sel_prob_map[row, col, image_idx] = beta

β 是一个隐藏马尔可夫模型 (HMM) 的后向消息，编码了当前像素下方的所有像素是否认为该源图是"好的匹配"。这为后续的源图采样提供了全局上下文。

后向消息暂时存储在 sel_prob_map_ 中（复用内存），在后续 sweep 逐行处理时被前向消息替换。

步骤 ②：空间传播

源码位置: patch_match_cuda.cu:1013

将上方像素 (row-1, col) 的 depth + normal 传播到当前行：

propagated_depth = PropagateDepth(
    prev_depth, prev_normal, row-1, row);

PropagateDepth 的计算：

上方像素 (row-1, col) 的 depth + normal 定义了一个 3D 平面:
  平面方程: n · (X - P) = 0
  其中 P = depth * 视线方向 (row-1 处的 3D 点)

当前像素 (row, col) 的视线 ray 与该平面求交:
  t = (n · P) / (n · ray_direction)
  propagated_depth = t * (ray_direction[2])  // 取 z 分量作为深度

关键：这不是简单复制邻域深度，而是通过平面几何正确传播斜面上的深度值。例如地面倾斜时，深度沿列方向的变化会被正确传播。

步骤 ③：随机扰动

源码位置: patch_match_cuda.cu:1021-1028

在当前最优假设基础上加随机扰动：

// 深度扰动: 在 [(1-p)*depth, (1+p)*depth] 均匀采样
rand_depth = PerturbDepth(perturbation, curr_depth, &rand_state);

// 法线扰动: 随机旋转矩阵，角度范围 [-p*π/2, p*π/2]
PerturbNormal(row, col, perturbation * M_PI, curr_normal, &rand_state, rand_normal);

扰动幅度 perturbation 按指数衰减：

iter=0, sweep=0:  p = 1.000    (100% 大范围搜索)
iter=0, sweep=1:  p = 0.841
iter=0, sweep=2:  p = 0.707
iter=0, sweep=3:  p = 0.595
iter=1, sweep=0:  p = 0.500
iter=2, sweep=0:  p = 0.250
iter=4, sweep=3:  p = 0.037    (3.7% 精细微调)

实现从粗到精的搜索策略：早期大范围探索解空间，后期精细调整。

深度扰动：乘性 vs 加性

深度扰动采用乘性范围 [(1-p)*depth, (1+p)*depth] 而非加性 [depth-δ, depth+δ]：

• 物理意义：深度测量的相对精度通常恒定（如 ±1%），而非绝对精度恒定（如 ±0.1m）。乘性扰动在不同深度范围上产生成比例的搜索范围
• 实际效果：depth=100m 时搜索 [50, 150]m（iter=0），depth=1m 时搜索 [0.5, 1.5]m。加性扰动 ±50m 在近距离时完全不合理

法线扰动：旋转矩阵

法线扰动使用旋转矩阵 R = Rx(α) × Ry(β) × Rz(γ)，角度 α, β, γ 在 [-perturbation/2, perturbation/2] 均匀采样。如果旋转后法线翻转（不再面向相机），则将扰动减半重试（最多 3 次）。

旋转矩阵扰动保证了法线始终在流形上（单位球面上），比直接在 (nx, ny, nz) 上加噪声后归一化更均匀。

指数衰减的数学保证

perturbation = 1/2^(iter + sweep/4)

• 与模拟退火的类比：perturbation 对应"温度"，指数衰减对应冷却调度
• 指数衰减 vs 线性衰减：线性衰减 p = 1 - step/total 在后期仍有较大残余扰动（如最后 10% 步骤的 p 从 0.1 线性降到 0）。指数衰减在最后几步的扰动极小（p < 0.04），保证收敛精度
• 收敛保证：perturbation 最终趋于 0 → 随机扰动退化 → 算法退化为纯坐标下降 → 保证收敛到局部最优

步骤 ④：源图选择概率计算

源码位置: patch_match_cuda.cu:1036-1072

为每张源图计算选择概率，由 4 个因子相乘：

因子 A：HMM 前向-后向消息

源码位置: patch_match_cuda.cu:664-798 (LikelihoodComputer)

HMM 状态空间 = {被选中, 未选中}，两个状态之间的转移概率：

自转移概率 = 0.99999  (极高的状态保持概率)
状态切换概率 = 0.00001

发射概率：

• 被选中状态: NCC 似然 = exp(-cost² / (2σ_ncc²)) × 归一化因子
• 未选中状态: 均匀概率 = 0.5

前向消息 α 在 sweep 过程中逐行更新，后向消息 β 已预计算。合并后：

sel_prob = prev_weight × prev_prob + (1 - prev_weight) × (α × β / total)

prev_weight 随迭代线性增长（0→0.95），给历史选择概率越来越大的权重，增加时序稳定性。

因子 B：三角化角概率

cos_angle = -dot(S, point) / (|S| × |point|)
S = source_camera_center - 3D_point

if cos_angle > cos(min_angle):   // 角度太小
    prob = 二次惩罚, 映射到 [0, 1]
else:
    prob = 1.0                   // 角度足够大，无惩罚

几何原理：三角化角是参考相机、3D 点、源相机构成的三角形的顶角。当该角度接近 0° 时，参考和源相机几乎从同一方向观察该点，深度估计的精度急剧下降：

深度不确定度 ∝ 特征定位误差 / sin(triangulation_angle)

  角度 30°: sin(30°) = 0.5   → 深度误差 = 2× 特征误差
  角度 5°:  sin(5°)  = 0.087 → 深度误差 ≈ 11.5× 特征误差
  角度 1°:  sin(1°)  = 0.017 → 深度误差 ≈ 57× 特征误差

二次惩罚而非硬阈值：硬阈值在边界处会导致源图突然切换，引入不连续性。二次惩罚在阈值附近平滑过渡，低于阈值后概率逐渐下降到 0。

因子 C：入射角概率

cos_incident = dot(source_ray_direction, normal)
prob = exp(-(1 - max(0, cos_incident))² / (2σ_incident²))

Foreshortening 效应：当源相机以大角度（接掠射）观察表面时：

• 表面上的小区域在图像中被压缩到更少的像素（foreshortening），匹配窗口覆盖的实际面积减小
• 纹理分辨率降低，特征区分性下降
• 非朗伯特表面（glossy、specular）的反射特性随角度变化更大

高斯衰减函数 exp(-(1-cosθ)²/(2σ²)) 在 θ=0°（正对表面）时为 1.0，θ=45° 时约 0.78（σ=0.9），θ=90°（掠射）时趋近 0。

因子 D：分辨率先验

计算单应变换后源图 patch 与参考图 patch 的面积比：

if 源图 patch 更小 (zoomed in):
    ratio = src_area / ref_area
else:
    ratio = ref_area / src_area

分辨率匹配的重要性：homography 将参考图的 NCC 窗口 warp 到源图后，patch 大小可能改变：

• 源图离场景更近 → 源图 patch 更大 → 分辨率更高 → 匹配更精细（好）
• 源图离场景更远 → 源图 patch 更小 → 分辨率不足 → 细节丢失（差）

取 min(src/ref, ref/src) 而非 src/ref 是因为两个方向的降质都应惩罚：无论是源图分辨率不足还是过度放大（可能引入插值伪影），匹配质量都会下降。比值接近 1.0 表示两视角的分辨率匹配良好。

合并:

sampling_prob[i] = sel_prob[i] × tri_prob[i] × inc_prob[i] × res_prob[i]

然后转换为 CDF 用于加权采样。

步骤 ⑤：Monte Carlo 源图采样 + 5 假设评估

源码位置: patch_match_cuda.cu:1094-1139

同时评估 5 个假设：

假设编号	depth 来源	normal 来源	说明
0	当前最优	当前最优	复用已存储 photo cost；几何模式仍需计算 geom cost
1	传播（上方）	传播（上方）	空间传播候选
2	随机扰动	随机扰动	探索新区域
3	当前最优	随机扰动	在已知好深度处探索新法线
4	随机扰动	当前最优	在已知好法线处探索新深度

采样过程（重复 num_samples=15 次）：

for sample = 0 to num_samples-1:
    # 按选择概率加权采样一张源图
    rand = random()
    src_idx = 第一个 CDF 超过 rand 的源图

    # 假设 0: 复用存储的 cost
    costs[0] += cost_map[row, col, src_idx]

    # 假设 1-4: 重新计算 NCC 代价
    for hypo = 1 to 4:
        costs[hypo] += ComputeNCC(depths[hypo], normals[hypo], src_idx)

采样结束后，累计代价最小的假设胜出：

best = argmin(costs[0..4])
depth_map[row, col]  = depths[best]
normal_map[row, col] = normals[best]

5 假设策略的设计原理

从坐标下降 (Coordinate Descent) 的视角理解：同时优化 depth 和 normal 构成一个 4 维搜索问题（1D depth + 3D normal），直接搜索代价高昂。5 个假设实现了对这个 4D 空间的高效覆盖：

策略	假设	优化维度	目的
保持	0: (curr_d, curr_n)	无	exploitation：保留已知的局部最优
传播	1: (prev_d, prev_n)	全维度	利用邻域的空间连贯性
探索	2: (rand_d, rand_n)	全维度	exploration：跳出局部最优
分离优化 depth	3: (rand_d, curr_n)	仅 depth	depth 单独优化（normal 已收敛时）
分离优化 normal	4: (curr_d, rand_n)	仅 normal	normal 单独优化（depth 已收敛时）

假设 3+4 的必要性：考虑一个 depth 已收敛到正确值但 normal 尚未收敛的像素。假设 0/1/3/4 都保持正确 depth，只有假设 2 和 4 探索新 normal。如果没有假设 4（只有 0-3），normal 的优化只能通过假设 2（全随机）进行，效率很低。分离优化允许在一个维度已收敛时高效优化另一个维度。

Monte Carlo 采样的自适应性

Monte Carlo 采样 num_samples=15 次（而非评估所有 N_src=20 张源图）的关键优势是自适应性：

• 概率集中时（如某源图 sel_prob = 0.8）：该源图被采样的期望次数 = 15 × 0.8 = 12 次，等价于主要用该最优源图评估 → 结果接近最优
• 概率均匀时（如 5 张源图各 0.2）：每张被采样的期望次数 = 15 × 0.2 = 3 次 → 均匀探索
• 对比 Top-K 选择：固定取概率最高的 K 张源图，无法适应不同分布。当 K 张源图概率接近时浪费计算，当某张源图概率远高于其他时又不够集中

采样概率 → CDF → 随机数 ∈ [0,1) → 第一个 CDF 超过随机数的源图

这种逆变换采样 (Inverse Transform Sampling) 保证了采样频率严格正比于选择概率。

步骤 ⑥：更新前向消息

源码位置: patch_match_cuda.cu:1154-1173

选定最优假设后，用新的 depth/normal 重新计算每张源图的 NCC 代价，更新前向消息 α 和 sel_prob_map。

步骤 ⑦：过滤（最后一轮 sweep）

源码位置: patch_match_cuda.cu:1175-1242

仅在最后一次迭代的最后一个 sweep 启用。统计每张源图的一致性：

for each source image:
    if NCC cost < filter_min_ncc:
        consistent_count++

if consistent_count < filter_min_num_consistent:
    depth_map[row, col] = 0    // 标记为无效

几何一致性模式下，额外检查前向-后向重投影误差。

光度一致性代价：NCC 计算细节

源码位置: patch_match_cuda.cu:460-558 (PhotoConsistencyCostComputer)

单应变换 (Homography)

给定像素 (row, col) 的 depth d 和 normal n，计算从参考图到源图的平面单应矩阵：

H = K_src × (R - T · n' / d) × K_ref⁻¹

其中:
  d = depth × (n · 视线方向)     // 平面到相机原点的距离
  R, T = 参考图到源图的相对旋转和平移
  n' = n 在相机坐标系下的表示

为什么 depth + normal 能定义 3D 平面

在相机坐标系下，像素 (row, col) 的视线方向为：

v = (inv_fx × col - cx/fx, inv_fy × row - cy/fy, 1)ᵀ

沿该视线方向、深度为 depth 的 3D 点为：

P = depth × v = (X, Y, depth)ᵀ

该点处的法线 n = (nx, ny, nz) 定义了一个过 P 的平面：

平面方程: nᵀ · (X - P) = 0
即: nx·X + ny·Y + nz·Z = nᵀ · P

平面上所有点满足 nᵀ · Q = d，其中 d = nᵀ · P = depth × (n · v) 就是平面到相机原点的有符号距离。由于法线面向相机 (n · v < 0)，d < 0，表示平面在相机前方。

平面诱导 Homography 的推导

两个相机之间的几何关系由旋转 R 和平移 T 描述。对于参考图像素 p_ref，其对应的 3D 点为 Q = K_ref⁻¹ × p_ref × Z（Z 是该点的深度）。如果 Q 落在平面 nᵀQ = d 上，则：

Z = d / (nᵀ × K_ref⁻¹ × p_ref)

代入投影方程 p_src ∝ K_src × (R × Q + T) 并展开：

p_src ∝ K_src × (R × K_ref⁻¹ × p_ref × d/(nᵀK_ref⁻¹p_ref) + T)
       ∝ K_src × (R - T × nᵀ/d) × K_ref⁻¹ × p_ref
       = H × p_ref

因此 平面上的所有像素对 (p_ref, p_src) 由一个 3×3 矩阵 H 关联，这是 homography 的定义。关键洞察：一般的双视图几何需要基础矩阵 F（非平面场景），但当我们假设局部平面时，关系退化为精确的 homography，不再需要逐像素搜索深度，只需指定平面参数 (depth, normal) 即可确定所有对应关系。

增量计算

warped 坐标通过增量计算避免重复矩阵乘法：

tform_base = H × [col, row, 1]ᵀ
tform_step_row = H × [0, 1, 0]ᵀ
tform_step_col = H × [1, 0, 0]ᵀ

逐像素遍历 NCC 窗口时，坐标递推为：

warped(row+1, col) = warped(row, col) + tform_step_row
warped(row, col+1) = warped(row, col) + tform_step_col

数值稳定性：图像坐标在千量级（如 2048），远大于颜色值 [0,1]。从 base 位置增量累加比每像素独立计算 H × [c,r,1] 更稳定，因为避免了 大数 × 小数 + 大数 × 大数 的精度损失。每个 NCC 窗口只需一次完整矩阵向量乘法 + 最多 window_size² 次加法。

共享内存参考图像

源码位置: patch_match_cuda.cu:343-414 (LocalRefImage)

每个线程块在共享内存中缓存一列参考图像 patch：

tile 尺寸: kNumRows × (3 × THREADS_PER_BLOCK) floats

第一行: 加载完整 tile (所有线程并行)
后续行: 上移一行，只加载最底部新行

NCC 公式

归一化互相关 (Normalized Cross-Correlation) 衡量两个图像 patch 的线性相关程度：

          Σ (I_ref - μ_ref) × (I_src - μ_src)
NCC = ─────────────────────────────────────────
       sqrt(Σ(I_ref - μ_ref)²) × sqrt(Σ(I_src - μ_src)²)

等价形式（编程更方便）：

NCC = Covar(ref, src) / sqrt(Var(ref) × Var(src))

数值含义：

NCC 值	含义
+1	完美正相关——两个 patch 亮度变化模式完全一致
0	无关——两个 patch 毫无关联
-1	完美负相关——亮度变化模式完全相反

NCC 通过减均值、除标准差消除了亮度和对比度差异。假设源图比参考图亮 50、对比度放大 2 倍（I_src = 50 + 2·I_ref），则 (I_src - μ_src) 依然是 (I_ref - μ_ref) 的 2 倍，但除以标准差后完全抵消，NCC = 1。

COLMAP 使用 cost = 1 - NCC，范围 [0, 2]，0 = 完美匹配，2 = 最差。

为什么选择 NCC 而非 SSD/SAD

不同视角的图像可能存在仿射亮度变化 I_src = a + b × I_ref（由曝光差异、传感器增益、gamma 校正等引起）。

• SSD（Sum of Squared Differences）对这种变化敏感：SSD = Σ(a + b·I₁ - I₂)² 随 a, b 变化而改变
• NCC 具有仿射不变性：

若 I_src = a + b × I_ref，则:
  E[I_src] = a + b × E[I_ref]
  Var(I_src) = b² × Var(I_ref)
  Covar(I_ref, I_src) = b × Var(I_ref)
  NCC = b × Var(ref) / sqrt(Var(ref) × b² × Var(ref)) = sign(b)

当 b > 0（大多数情况）时，NCC = 1（完美相关），与 a, b 的具体值无关。这使得 NCC 在不同曝光条件下保持稳定的匹配质量。

双边权重的物理意义

普通加权求和给每个像素一个权重 Σ w(i) × I(i)。双边就是从两个维度同时决定权重：

w(i) = exp(-距离²/(2σ_s²)) × exp(-颜色差²/(2σ_c²))
       ─────────────────────   ────────────────────────
          空间权重                颜色权重

• 空间权重 exp(-dist²/(2σ_s²))：高斯衰减，σ_s = window_radius 时窗口边缘权重约 0.14。给中心像素更高权重，同时避免窗口边缘像素主导统计量。

• 颜色权重 exp(-Δcolor²/(2σ_c²))：σ_c = 0.2（灰度范围 [0,1]）。关键作用是跨遮挡边界的自适应降权：

  • 当窗口跨越遮挡边界时，边界另一侧的像素颜色差异大 → 权重低 → 不参与统计
  • 这避免了遮挡物像素"污染"NCC 计算，使得遮挡边界处的深度估计更锐利
  • 等效于只在颜色相似（大概率属于同一平面）的像素上计算 NCC

遮挡边界示例：

场景: 一面白墙前面有根深色柱子，NCC 窗口中心在墙上（白色）

┌─────────────────┐
│ 白 白 白 黑 黑   │  ← 窗口内一半是墙，一半是柱子
│ 白 [白] 白 黑 黑 │  ← 中心是墙
│ 白 白 白 黑 黑   │
└─────────────────┘

普通高斯加权: 黑色像素距离近 → 权重高 → 混入柱子颜色 → NCC 被拉低
双边加权:     黑色像素和中心(白色)颜色差大 → 权重趋近0 → 只用白色像素算 NCC

与普通高斯滤波的对比：

	高斯滤波	双边滤波
空间距离	远了降权 ✓	远了降权 ✓
颜色差异	不管	差了降权 ✓
跨边界行为	会模糊 ✗	保持锐利 ✓

预计算的数学依据

参考图像侧的统计量在迭代中不变，因为参考图像本身不随 depth/normal 变化而改变：

ref_sum = Σ w(i) × ref_color(i)          ← 常量（参考图不变）
ref_sq_sum = Σ w(i) × ref_color(i)²      ← 常量

只有源图像侧的 warp 坐标（以及读出的颜色值）随假设变化。因此 sum_image 和 squared_sum_image 只需在初始化时计算一次，每次假设评估只更新源图侧的 src_color_sum、src_color_sq_sum、src_ref_color_sum。这将 NCC 的每像素计算量减少了约一半。

具体计算

ref_color_var = squared_sum_image - sum_image²        (预计算)
src_color_sum = Σ w(i) × src_color(i)
src_color_sq_sum = Σ w(i) × src_color(i)²
cross_sum = Σ w(i) × ref_color(i) × src_color(i)

src_color_var = src_color_sq_sum / W - (src_color_sum / W)²
covar = cross_sum / W - (ref_sum / W) × (src_color_sum / W)

NCC = covar / sqrt(max(ref_color_var, ε) × max(src_color_var, ε))

其中 W = Σ w(i) 是权重和，ε = 1e-5 防止除零。若任一方差 < ε，直接设 cost = 2.0（最差代价）。

代价范围: [0, 2.0]，其中 0 = 完美匹配，2.0 = 最差。

HMM 源图选择的数学模型

源码位置: patch_match_cuda.cu:664-798 (LikelihoodComputer)

设计动机

MVS 重建中，并非所有源图对每个像素都有用。一张源图可能因为遮挡、视角过大、分辨率不匹配等原因对某些像素无效。暴力评估所有源图既浪费计算又可能引入坏匹配。

核心观察：一张源图在相邻像素处的"好坏"状态通常是空间连续的。例如，源图 A 在参考图像的某个区域可见（匹配好），在另一个区域被物体遮挡（匹配差）。遮挡边界处状态突变，但边界两侧各自连续。

HMM 建模思路：将 sweep 方向（如从上到下）的像素序列视为时间轴，每张源图的"好坏"状态形成一条沿空间轴的一维 HMM。前向-后向消息传递让每个像素的源图选择不仅基于本地 NCC，还融合了整列像素的全局证据。

遮挡场景示例 (sweep 方向: 上→下):

  行 0-50:   源图 A 可见，NCC 高 → sel_prob 高 → 被频繁采样
  行 51:     遮挡边界，NCC 骤降 → 前向消息开始下降
  行 52-100: 源图 A 被遮挡，NCC 低 → sel_prob 低 → 很少被采样

  HMM 的作用: 行 51 处 sel_prob 不是瞬间跳变，而是平滑过渡
             （因为 0.99999 的自转移概率使状态变化需要多行证据累积）

模型定义

对每个像素，每张源图有一个隐藏状态 s ∈ {0, 1}：

• s = 1: 该源图被选中（匹配质量好）
• s = 0: 该源图未被选中（匹配质量差或不一致）

转移模型

P(s_t = s_{t-1}) = 0.99999    (极高的状态保持概率)
P(s_t ≠ s_{t-1}) = 0.00001    (极低的状态切换概率)

其中 t 是行索引（在垂直 sweep 中）。

发射模型

P(cost | s=1) = exp(-cost² / (2σ_ncc²)) × Z      (好的匹配 → 高 NCC → 低 cost → 高概率)
P(cost | s=0) = 0.5                                (均匀分布)

归一化因子 Z = 2 / (√(2π)·σ·erf(√2/σ))，确保发射概率在 cost∈[0,2] 上合理归一化。其中 cost = 1 - NCC，NCC∈[-1,1]，故 cost∈[0,2]。

前向-后向消息传递

前向消息 α 的直觉：在从上到下扫描时，α 编码了"从第 0 行到当前行，所有像素的证据是否支持选中该源图"。它像一个从上往下传递的"信任分数"——如果上方很多像素都认为该源图好，α 会维持在高位；如果遇到连续多行低 NCC，α 逐渐下降。

前向消息 α (sweep 方向内逐行更新):

α_t(1) ∝ P(cost_t | s=1) × [0.99999 × α_{t-1}(1) + 0.00001 × α_{t-1}(0)]
α_t(0) ∝ P(cost_t | s=0) × [0.00001 × α_{t-1}(1) + 0.99999 × α_{t-1}(0)]

后向消息 β 的直觉：β 编码了"从当前行到底部，该源图是否被看好"。在 sweep 开始前从底部向上计算，β 让当前像素在决策时不仅看"上方说了什么"（α），还知道"下方会说什么"（β）。

为什么用 0.99999 而非 0.9：转移概率 P(保持) = 0.99999 意味着：

• 状态切换需要约 1/0.00001 = 100,000 个像素的证据才能翻转
• 在实际图像（高度通常 < 4000 像素）中，这意味着一次遮挡不会导致选择概率剧烈波动
• 只有持续多行的强证据（连续高/低 NCC）才会改变状态
• 效果：选择概率具有空间平滑性，避免单像素噪声导致源图选择抖动

后向消息 β (逆 sweep 方向预计算):

β_t(1) ∝ Σ_s' P(s_{t+1}=s' | s_t=1) × P(cost_{t+1} | s') × β_{t+1}(s')
β_t(0) ∝ Σ_s' P(s_{t+1}=s' | s_t=0) × P(cost_{t+1} | s') × β_{t+1}(s')

合并:

P(s=1 | all observations) = α_t(1) × β_t(1) / [(1-α_t(1))×(1-β_t(1)) + α_t(1)×β_t(1)]

其中 α_t 是标量（前向消息 P(s=1|o_0…o_t)），β_t 也是标量（后向消息 P(s=1|o_t…o_{H-1})）。代码中 ComputeSelProb 函数直接实现了这个公式。

时序平滑

相邻 sweep 方向之间的选择概率通过线性插值平滑：

final_prob = prev_weight × prev_prob + (1 - prev_weight) × (α × β / total)
prev_weight = (iter×4 + sweep) / (num_iterations×4)

跨迭代的信息累积：每个 sweep 方向计算的选择概率 α×β/(total) 可能因当轮 depth/normal 的随机性而抖动。通过将上一轮的结果 prev_prob 以递增权重混合进来，算法在迭代过程中逐步建立起稳定的源图选择共识。

prev_weight 从 0 线性增长到 (N_iter×4 - 1) / (N_iter×4)（默认 5 次迭代时最大为 19/20 = 0.95）：

• 早期 (prev_weight ≈ 0)：完全信任当前 sweep 的计算，允许快速适应新的 depth/normal
• 后期 (prev_weight ≈ 0.95)：主要信任历史累积概率，depth/normal 已基本收敛，不需要源图选择剧烈变化
• 始终保留 5% 新信息：即使到最后也不完全锁定，保留对局部异常的响应能力

深度与法线参数化

为什么用 (depth, normal) 而非其他参数化

MVS 算法需要对每个像素估计一个 3D 位置。常见的参数化选择：

参数化	自由度	优势	劣势
Disparity (视差)	1	正面平面场景下线性化	无法表示倾斜平面，需要视角已知
Depth (深度)	1	简单直观	只能表示正面朝向的平面，倾斜面需要逐像素独立深度
3D 坐标 (X,Y,Z)	3	最通用	参数空间不平滑，优化困难，需要更多样本
Depth + Normal	4	局部平面假设，参数空间平滑	比纯深度多 3 个参数

(depth, normal) 参数化的核心优势是局部平面假设：

1. 真实世界的表面在局部范围内近似平面（墙壁、地面、桌面等）
2. depth + normal 定义一个过该像素 3D 点的平面，邻域像素共享同一平面假设
3. 空间传播时，邻域的平面假设可以直接通过几何传播（而非启发式插值）
4. 相比逐像素独立深度，平面假设通过法线提供了深度在空间上的一阶梯度信息，大幅减少需要的随机搜索次数

深度

每像素一个 float，表示沿光轴的深度值。3D 点重建：

X = depth × (inv_fx × col - cx/fx)
Y = depth × (inv_fy × row - cy/fy)
Z = depth

法线

每像素一个 3D 单位向量 (nx, ny, nz)，约束面向相机：

dot(normal, view_ray_direction) < 0

depth + normal 共同定义该像素处的局部 3D 平面。

平面传播

空间传播时，上方像素的 (depth, normal) 定义了一个 3D 平面。实现上 PropagateDepth 函数做了 2D 简化：只在 (row方向, depth) 平面内做线段求交，仅使用法线的 ny 和 nz 分量（忽略 nx），因为 sweep 方向沿列传播，横向分量不影响交点深度。

2D 简化（row-depth 平面）:
  上方像素的 ray 终点: P1 = (depth1 * (inv_fy*row1 - cy/fy), depth1)
  平面法线方向:        dir = (nz, -ny)

  当前像素的 ray 方向: (inv_fy*row2 - cy/fy, 1)

  线段 ((0,0)→ray2方向) 与 线段 (P1→P1+dir) 的交点:
    t = (y1*x2 - x1*y2) / ((x2-x1) + ray2_x * (y1-y2))

  propagated_depth = t

概念上等价于完整 3D ray-plane 求交，但避免了 3D 运算和 nx 分量，计算更高效。

迭代调度

源码位置: patch_match_cuda.cu:1359-1512

完整时间表

Step 0: ComputeInitialCost
        用初始随机 depth/normal 计算所有源图的 NCC 代价

Step 1-20: 5 iterations × 4 sweeps
  ├─ iter 0, sweep 0: p=1.000, prev_w=0.00  (大范围搜索)
  ├─ iter 0, sweep 1: p=0.841, prev_w=0.05
  ├─ iter 0, sweep 2: p=0.707, prev_w=0.10
  ├─ iter 0, sweep 3: p=0.595, prev_w=0.15
  ├─ iter 1, sweep 0: p=0.500, prev_w=0.20
  ├─ iter 1, sweep 1: p=0.420, prev_w=0.25
  ├─ ...
  ├─ iter 4, sweep 2: p=0.035, prev_w=0.90
  └─ iter 4, sweep 3: p=0.037, prev_w=0.95  (精细微调 + 过滤)

扰动衰减曲线

perturbation = 1.0 / 2^(iter + sweep/4)

迭代	Sweep 0	Sweep 1	Sweep 2	Sweep 3
0	100.0%	84.1%	70.7%	59.5%
1	50.0%	42.0%	35.4%	29.7%
2	25.0%	21.0%	17.7%	14.9%
3	12.5%	10.5%	8.8%	7.4%
4	6.3%	5.3%	4.4%	3.7%

旋转策略

源码位置: patch_match_cuda.cu:1810

每次 sweep 后，整个 GPU 状态旋转 90°：

旋转的数据:
  - depth_map, normal_map, cost_map
  - sel_prob_map, prev_sel_prob_map
  - ref_image, sum_image, squared_sum_image
  - rand_state_map (随机数状态)

切换的资源:
  - poses_texture_[rotation]  (4 种预计算的旋转位姿)
  - ref_K[rotation], ref_inv_K[rotation] (CUDA 常量内存)

4 个 sweep 的方向：

Sweep	原始方向	旋转后实际遍历
0	上→下	上→下
1	上→下 (旋转 90°)	右→左
2	上→下 (旋转 180°)	下→上
3	上→下 (旋转 270°)	左→右

设计理由：通过旋转数据结构而非改变遍历方向，复用完全相同的 CUDA kernel，避免为不同方向写不同代码。

几何一致性阶段 (Phase 2)

两阶段执行架构

源码位置: patch_match.cc:170-207 (PatchMatchController::Run())

PatchMatchController 采用两阶段执行策略：

if geom_consistency == true:
    Pass 1 (Photometric):
        options.geom_consistency = false
        options.filter = false
        → 对所有图像运行纯光度 PatchMatch
        → 输出标记为 "photometric" 类型
        → depth/normal 随机初始化

    Pass 2 (Geometric):
        options = 原始选项 (geom_consistency=true, filter=true)
        input_type = "photometric" (读取 Pass 1 的输出)
        → 对所有图像运行几何一致性 PatchMatch
        → 输出标记为 "geometric" 类型
        → depth/normal 从 photometric 结果初始化（非随机）

else:
    单阶段执行，filter 按配置决定

为什么需要两遍——源图的深度从哪来：

重投影误差的计算要求源图也有深度：参考图用自己深度反投影到 3D → 投影到源图 → 要用源图的深度反投影回来。但源图的深度从哪来？

答案就是 Pass 1 已经对所有图像都算了一遍深度图：

Pass 1 结束后，每张图都有自己的深度图:

  depth_maps/
    image_001.photometric.bin   ← 图像 1 的深度图
    image_002.photometric.bin   ← 图像 2 的深度图
    image_003.photometric.bin   ← 图像 3 的深度图
    ...

重投影时:

  参考图 (500,300)  depth=10m
      ↓ 反投影
    3D 点 A
      ↓ 投影
  源图 (480,290)    ← 落到源图某个像素
      ↓ 读这个像素的深度图值
  源图 (480,290) 的 depth = 9.8m  ← Pass 1 算出来的
      ↓ 用源图自己的深度反投影
    3D 点 A'
      ↓ 投影回参考图
  参考图 (498, 295)  ← 偏了 2.5 像素

Pass 1 让每张图都先猜好深度，Pass 2 让它们互相用对方的深度来验证自己猜得对不对。

关键设计：Pass 1 不启用过滤（filter=false），确保 photometric 结果尽可能完整；Pass 2 的 depth_min/depth_max 从 sparse SfM 自动计算；sigma_spatial 如果 ≤ 0，自动设为 window_radius。

自适应过滤阈值

源码位置: patch_match.cc:458-460

控制器会自适应调整 filter_min_num_consistent：

patch_match_options.filter_min_num_consistent =
    std::min(static_cast<int>(used_image_idxs.size()) - 1,
             patch_match_options.filter_min_num_consistent);

当某问题可用的源图数量少于 filter_min_num_consistent + 1 时，自动降低阈值以避免过滤掉所有像素。

初始化

depth 和 normal 直接使用 Photometric 阶段的结果（非随机初始化）。

额外代价项

在 NCC 代价外加入前向-后向重投影误差：

geom_cost = Σ_src  (forward_reproj_error + backward_reproj_error)

总代价 = photo_cost + λ × geom_cost
其中 λ = geom_consistency_regularizer = 0.3

几何代价项直接以权重 λ 加到 光度代价上，而非 (1-λ)×photo + λ×geom 的加权平均。photo_cost 范围 [0, 2]，geom_cost 截断在 max_cost (3.0px)，故 λ×geom_cost 最大贡献 0.3×3=0.9。

重投影误差的直觉：去了一趟能不能回到原点

两张照片拍同一面墙，墙上有个标记点 A。重投影误差就是：从参考图出发，走 参考图→3D→源图→3D→参考图 一圈，看看能不能回到出发点。

深度正确时——回到原点：

参考图像素 (500,300)，真实深度 10m

  参考图 (500,300)                    源图 (480,290)
       │                                  ▲
       │ ①反投影: depth=10m               │
       ▼                                  │
     3D 点 A（在墙上）──②投影到源图──→│    ← 恰好是 A 在源图的位置
                                          │
       ▲                                  │ ③读源图深度=10m
       │                                  │ ④反投影: 还是 A
       │           ⑤投影回参考图─────────┘
  回到 (500,300) = 出发点

  误差 = 0 像素 ✓

深度猜错时——回不到原点：

参考图像素 (500,300)，猜深度 5m（真实是 10m）

  参考图 (500,300)                    源图 (460,275)
       │                                  ▲
       │ ①反投影: depth=5m               │
       ▼                                  │
     3D 点在墙前面 ──②投影到源图──→│     ← 不是 A，跑到别处去了！
  （空中，不在墙上）                        │
                                          │ ③读这个位置的深度=8m
       ▲                                  │ ④反投影: 墙上另一个点
       │           ⑤投影回参考图─────────┘
  回到 (470,280) ≠ 出发点 (500,300)

  误差 = √((500-470)² + (300-280)²) = √(900+400) ≈ 36 像素 ✗

核心逻辑：

正确的深度:  参考图 → 3D → 源图 → 3D → 参考图 = 回到原点   误差=0
错误的深度:  参考图 → 3D → 源图 → 3D → 参考图 ≠ 回到原点   误差>0

两张图互相验证：我用我的深度算出 3D 点去你那边看，你用你的深度从那边算回来。如果对不上（回不到原点），说明至少有一方的深度有问题。

前向投影的具体计算步骤

步骤 1: 反投影（参考图像素 → 3D 点）
  3D_point = depth × (inv_fx×col - cx/fx, inv_fy×row - cy/fy, 1)

步骤 2: 投影到源图（3D 点 → 源图像素，使用投影矩阵 P）
  src_col = (P[0]*X + P[1]*Y + P[2]*Z + P[3]) / (P[8]*X + P[9]*Y + P[10]*Z + P[11])
  src_row = (P[4]*X + P[5]*Y + P[6]*Z + P[7]) / (P[8]*X + P[9]*Y + P[10]*Z + P[11])

步骤 3: 读取源图该位置的深度
  src_depth = src_depth_map[src_row, src_col]

步骤 4: 反投影回 3D（源图像素 → 3D 点，使用逆投影矩阵 inv_P）
  backward_point = inv_P × (src_col × src_depth, src_row × src_depth, src_depth)

步骤 5: 投影回参考图（3D 点 → 参考图像素）
  backward_col = fx × backward_X / backward_Z + cx
  backward_row = fy × backward_Y / backward_Z + cy

步骤 6: 计算重投影误差
  error = min(max_cost, sqrt((col - backward_col)² + (row - backward_row)²))

当两张图的深度都正确时，round-trip 经过正确的 3D 点，返回原始像素坐标，误差为 0。当任一图的深度有误时，3D 点不匹配，投影回参考图后产生偏移，且偏移量与深度误差成正比。

λ 权重分析

photo_cost 范围:  [0, 2.0]   (1-NCC)
λ × geom_cost:   [0, 0.9]   (0.3 × max 3.0px)

几何项贡献占比:
  当 photo=0.1, geom=0.5px: total=0.1+0.15=0.25, 几何占 60%
  当 photo=0.5, geom=0.5px: total=0.5+0.15=0.65, 几何占 23%
  当 photo=1.0, geom=3.0px: total=1.0+0.9=1.9,   几何占 47%

λ=0.3 使几何项在光度匹配较好（photo < 0.3）时具有显著影响力，能纠正"NCC 看起来好但几何不一致"的情况（如重复纹理区域），但在光度匹配差时不会主导代价（避免误杀）。

过滤

最后一轮 sweep 使用更严格的过滤条件：

对每张源图:
  前置条件: 三角化角 > filter_min_triangulation_angle (3.0°)
           AND 入射角 > 0 (源图在该法线方向可见)

  光度一致性: sel_prob >= ComputeNCCProb(1 - filter_min_ncc)
             (HMM 后验选择概率超过阈值)
  AND
  几何一致性: 前向-后向重投影误差 < filter_geom_consistency_max_cost (1.0 px)

if 一致源图数 < filter_min_num_consistent:
    depth = 0  (标记为无效)

光度一致性过滤不是直接比较 NCC 代价值，而是比较 HMM 后验选择概率 sel_prob 与阈值 ComputeNCCProb(1 - filter_min_ncc)。这意味着过滤利用了整个 sweep 方向上的全局上下文信息，而非单像素代价。

性能特征

计算复杂度

阶段	复杂度	瓶颈
双边预过滤	O(W × H × window²)	参考图像处理
单次 Sweep	O(W × H × N_samples × NCC_window²)	NCC 计算
全部迭代	O(N_iter × 4 × W × H × N_samples × NCC_window²)	-

GPU 利用

• 纹理缓存：源图像存储在 CUDA layered texture 中，硬件双线性插值
• 共享内存：参考图像 patch 缓存在共享内存，减少全局内存访问
• 每列一线程：纵向扫描充分利用空间局部性和共享内存 tile

内存占用

以 2048×1536 图像、20 张源图为例：

数据	尺寸	内存
depth_map	2048×1536 × 4B	12 MB
normal_map	2048×1536 × 3 × 4B	36 MB
cost_map	2048×1536 × 20 × 4B	240 MB
sel_prob_map	2048×1536 × 20 × 4B	240 MB
src_images	20 × ~3MB	~60 MB
合计		~590 MB

关键设计总结

设计决策	方法	优势
光度代价	双边加权 NCC	对遮挡边界和颜色跳变鲁棒
空间传播	平面传播 (非像素复制)	正确处理斜面深度传递
源图选择	HMM 前向-后向消息	自适应选择最优视角，避免遮挡干扰
源图采样	Monte Carlo (15/20)	O(N_samples) 而非 O(N_src)
假设评估	5 候选并行	兼顾利用与探索
扰动策略	指数衰减	粗到精收敛
方向遍历	数据结构旋转	复用同一 kernel，代码简洁
几何约束	前向-后向重投影	过滤光度一致的几何异常值

参考资料

• Schönberger et al., “Pixelwise View Selection for Unstructured Multi-View Stereo”, ECCV 2016
• Bailer et al., “PatchMatch Stereo”
• COLMAP 源码: https://github.com/colmap/colmap
• https://colmap.github.io/

文章链接：
https://www.zywvvd.com/notes/3d/colmap-patchmatch-stereo/colmap-patchmatch-stereo/