1
2

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

1
2
3
4
5

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

1
2
3
4
5
6
7
8
9
10
11

COLMAP MVS 流水线:

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

1
2
3
4
5
6

[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)

1
2
3
4
5
6
7
8
9
10
11
12
13

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

1
2
3
4
5
6
7
8
9
10

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)   ← 同上

1
2
3
4
5
6
7
8

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°

1
2
3
4
5
6
7
8
9
10
11
12
13
14

第 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（差）

1
2
3
4
5
6
7
8
9
10
11
12

真实情况: 像素 (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 ✗

1
2
3
4
5

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

1
2
3
4
5
6
7
8

每轮 sweep 做两件事:

1. 传播: 抄邻居的好结果 → 正确信息向前推进一步
2. 随机扰动: 每个像素自己也在探索新深度 → 可能直接猜中正确答案
两个机制叠加:
  传播让正确值稳步扩散（确定性推进）
  随机扰动让远处像素有机会自己找到正确答案（不需要等传播传过来）

1
2
3
4
5

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

1
2
3
4

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

1
2
3
4
5
6
7
8

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

1
2
3
4
5
6

对每张源图 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

1
2

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

1
2
3
4
5
6
7

上方像素 (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 分量作为深度

1
2
3
4
5

// 深度扰动: 在 [(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);

1
2
3
4
5
6
7

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% 精细微调)

1

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

1
2

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

1

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

1
2
3
4
5
6
7

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                   // 角度足够大，无惩罚

1
2
3
4
5

深度不确定度 ∝ 特征定位误差 / sin(triangulation_angle)
  角度 30°: sin(30°) = 0.5   → 深度误差 = 2× 特征误差
  角度 5°:  sin(5°)  = 0.087 → 深度误差 ≈ 11.5× 特征误差
  角度 1°:  sin(1°)  = 0.017 → 深度误差 ≈ 57× 特征误差

1
2

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

1
2
3
4

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

1

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

1
2
3
4
5
6
7
8
9
10
11

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)

1
2
3

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

1

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

1
2
3
4
5
6

for each source image:
    if NCC cost < filter_min_ncc:
        consistent_count++
if consistent_count < filter_min_num_consistent:
    depth_map[row, col] = 0    // 标记为无效

1
2
3
4
5
6

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

1

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

1

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

1
2

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

1

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

1
2
3

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

1
2
3

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

1
2

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

1
2
3
4

tile 尺寸: kNumRows × (3 × THREADS_PER_BLOCK) floats
第一行: 加载完整 tile (所有线程并行)
后续行: 上移一行，只加载最底部新行

1
2
3

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

1

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

1
2
3
4
5

若 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)

1
2
3

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

1
2
3
4
5
6
7
8
9
10

场景: 一面白墙前面有根深色柱子，NCC 窗口中心在墙上（白色）
┌─────────────────┐
│ 白 白 白 黑 黑   │  ← 窗口内一半是墙，一半是柱子
│ 白 [白] 白 黑 黑 │  ← 中心是墙
│ 白 白 白 黑 黑   │
└─────────────────┘
普通高斯加权: 黑色像素距离近 → 权重高 → 混入柱子颜色 → NCC 被拉低
双边加权:     黑色像素和中心(白色)颜色差大 → 权重趋近0 → 只用白色像素算 NCC

1
2

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

1
2
3
4
5
6
7
8
9

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, ε))

1
2
3
4
5
6
7
8

遮挡场景示例 (sweep 方向: 上→下):
  行 0-50:   源图 A 可见，NCC 高 → sel_prob 高 → 被频繁采样
  行 51:     遮挡边界，NCC 骤降 → 前向消息开始下降
  行 52-100: 源图 A 被遮挡，NCC 低 → sel_prob 低 → 很少被采样
  HMM 的作用: 行 51 处 sel_prob 不是瞬间跳变，而是平滑过渡
             （因为 0.99999 的自转移概率使状态变化需要多行证据累积）

1
2

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

1
2

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

1
2

α_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)]

1
2

β_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')

1

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

1
2

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

1
2
3

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

1

dot(normal, view_ray_direction) < 0

1
2
3
4
5
6
7
8
9
10

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

1
2
3
4
5
6
7
8
9
10
11
12
13

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  (精细微调 + 过滤)

1

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

1
2
3
4
5
6
7
8
9

旋转的数据:
  - 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 常量内存)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

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 按配置决定

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

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 像素

1
2
3

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

1
2
3
4

geom_cost = Σ_src  (forward_reproj_error + backward_reproj_error)
总代价 = photo_cost + λ × geom_cost
其中 λ = geom_consistency_regularizer = 0.3

1
2
3
4
5
6
7
8
9
10
11
12
13
14

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14

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

1
2

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

步骤 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)²))

1
2
3
4
5
6
7

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%

1
2
3
4
5
6
7
8
9
10
11

对每张源图:
  前置条件: 三角化角 > 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  (标记为无效)