用Python和NumPy手把手实现光度立体法：从多张照片重建物体3D表面细节

想象一下，你只需要用手机环绕物体拍摄几张照片，就能自动生成这个物体的3D表面细节——包括每个微小凹凸的法线方向和材质反光特性。这就是光度立体法(Photometric Stereo)的魅力所在。不同于昂贵的3D扫描设备，这种方法仅需普通相机和可控光源，就能实现惊人的表面细节重建效果。

本文将带你用Python和NumPy一步步实现这个计算机视觉中的经典算法。无论你是想为游戏开发快速生成材质贴图，还是需要为工业检测分析产品表面缺陷，掌握这项技术都能大幅提升你的工作效率。我们将完全从实践角度出发，避开复杂的数学推导，专注于代码实现和实际问题的解决。

1. 准备工作与环境搭建

在开始编码前，我们需要确保准备好所有必要的工具和数据。光度立体法的核心输入是一组在相同视角、不同光照条件下拍摄的物体照片，以及对应的光源方向信息。

1.1 安装必要的Python库

我们将使用以下Python库来实现算法：

pip install numpy opencv-python matplotlib scipy

• NumPy ：处理所有矩阵运算
• OpenCV ：读取和处理图像
• Matplotlib ：可视化结果
• SciPy ：可选，用于后期法线图平滑处理

1.2 准备输入图像数据

理想的光度立体法输入图像应满足以下条件：

• 相机固定，仅光源位置变化
• 物体位置和姿态保持不变
• 使用漫反射为主的物体（避免镜面高光）
• 最好在暗室环境中拍摄，减少环境光干扰

假设我们已经拍摄了8张照片（建议最少3张，越多结果越稳定），存储在 ./images/ 目录下，命名为 img1.jpg 到 img8.jpg 。

import cv2
import numpy as np
import os

# 读取所有图像
image_files = sorted([f for f in os.listdir('./images') if f.endswith('.jpg')])
images = [cv2.imread(f'./images/{f}', cv2.IMREAD_GRAYSCALE) for f in image_files]
images = np.array(images, dtype=np.float32) / 255.0  # 归一化到[0,1]

print(f"加载了{len(images)}张图像，每张尺寸为{images[0].shape}")

1.3 光源方向测量

准确的光源方向对结果至关重要。测量方法包括：

1. 镜面球法 ：在场景中放置一个镜面球体，通过高光位置计算光源方向
2. 已知物体法 ：使用已知几何形状的物体（如圆柱体）反推光源
3. 专业测光设备 ：使用光度探头直接测量

假设我们已经测得8个光源方向，存储为一个8x3的NumPy数组：

light_directions = np.array([
    [0.1, -0.2, 1],   # 光源1方向(x,y,z)
    [0.3, -0.1, 1],   # 光源2方向
    # ... 其他光源方向
    [-0.2, 0.3, 1]    # 光源8方向
])

# 归一化为单位向量
light_directions = light_directions / np.linalg.norm(light_directions, axis=1, keepdims=True)

注意：光源方向向量应为从物体指向光源的方向，且z分量通常为正（假设相机沿z轴正方向观察）

2. 核心算法实现

有了图像和光源数据，我们现在可以实现光度立体法的核心计算部分。算法主要分为三步：计算法线和反射率、后处理、重光照演示。

2.1 计算法线和反射率

根据Lambertian反射模型，图像亮度可以表示为：

I = ρ * (n · s)

其中：

• I：观察到的像素亮度
• ρ：表面反射率（albedo）
• n：表面法线（单位向量）
• s：光源方向（单位向量）

对于每个像素，我们有多个方程（对应多张图像），可以解出ρ和n。

def compute_normals_and_albedo(images, light_directions, mask=None):
    """
    计算法线贴图和反射率图
    :param images: [k,h,w] k张h x w的图像
    :param light_directions: [k,3] k个光源方向
    :param mask: [h,w] 可选，背景掩码
    :return: normals [h,w,3], albedo [h,w]
    """
    k, h, w = images.shape
    S = light_directions  # [k,3]
    I = images.reshape(k, -1)  # [k, h*w]
    
    if mask is None:
        mask = np.ones((h, w), dtype=bool)
    mask_flat = mask.flatten()  # [h*w]
    
    # 初始化输出
    normals = np.zeros((3, h*w), dtype=np.float32)
    albedo = np.zeros(h*w, dtype=np.float32)
    
    # 只处理mask内的像素
    valid_pixels = I[:, mask_flat]  # [k, num_valid]
    
    # 最小二乘解 ρn = (S^T S)^-1 S^T I
    rho_n = np.linalg.pinv(S.T @ S) @ S.T @ valid_pixels
    
    # 计算反射率 ρ = ||ρn||
    valid_albedo = np.linalg.norm(rho_n, axis=0)
    
    # 计算法线 n = ρn / ρ
    valid_normals = rho_n / (valid_albedo + 1e-10)
    
    # 填充结果
    normals[:, mask_flat] = valid_normals
    albedo[mask_flat] = valid_albedo
    
    # 调整形状
    normals = normals.T.reshape(h, w, 3)
    albedo = albedo.reshape(h, w)
    
    # 确保法线是单位向量
    normals = normals / (np.linalg.norm(normals, axis=2, keepdims=True) + 1e-10)
    
    return normals, albedo

2.2 法线图后处理

原始计算的法线图可能存在噪声，我们可以应用一些后处理技术来改善结果：

def process_normals(normals, albedo, mask):
    """
    法线图后处理
    :param normals: [h,w,3] 原始法线
    :param albedo: [h,w] 反射率
    :param mask: [h,w] 背景掩码
    :return: 处理后的法线
    """
    from scipy.ndimage import median_filter
    
    # 1. 中值滤波去除离群点
    filtered_normals = normals.copy()
    for i in range(3):
        filtered_normals[..., i] = median_filter(filtered_normals[..., i], size=3)
    
    # 2. 只在mask区域内处理
    filtered_normals[~mask] = normals[~mask]
    
    # 3. 重新归一化
    norms = np.linalg.norm(filtered_normals, axis=2)
    filtered_normals = filtered_normals / (norms[..., np.newaxis] + 1e-10)
    
    return filtered_normals

2.3 重光照演示

有了法线和反射率，我们可以模拟物体在任何光源方向下的外观：

def relight_scene(albedo, normals, light_direction):
    """
    重新照亮场景
    :param albedo: [h,w] 反射率图
    :param normals: [h,w,3] 法线图
    :param light_direction: [3,] 光源方向（单位向量）
    :return: [h,w] 重光照图像
    """
    # 计算每个像素的亮度
    shading = np.sum(normals * light_direction, axis=2)
    shading = np.clip(shading, 0, 1)  # 限制在[0,1]范围
    
    # 应用反射率
    relit = albedo * shading
    
    return relit

3. 结果可视化与分析

计算完成后，我们需要有效展示和评估结果。法线图通常以RGB颜色表示，其中R、G、B通道分别对应法线的x、y、z分量。

3.1 可视化法线图

import matplotlib.pyplot as plt

def visualize_normals(normals, mask=None):
    """
    可视化法线图
    :param normals: [h,w,3] 法线图
    :param mask: [h,w] 可选背景掩码
    """
    # 将法线从[-1,1]映射到[0,1]用于显示
    display_normals = (normals + 1) / 2
    
    if mask is not None:
        # 将背景设为黑色
        display_normals[~mask] = 0
    
    plt.imshow(display_normals)
    plt.title('Normal Map')
    plt.axis('off')
    plt.show()

3.2 可视化反射率图

def visualize_albedo(albedo, mask=None):
    """
    可视化反射率图
    :param albedo: [h,w] 反射率
    :param mask: [h,w] 可选背景掩码
    """
    plt.imshow(albedo, cmap='gray')
    plt.title('Albedo Map')
    plt.axis('off')
    plt.colorbar()
    plt.show()

3.3 重光照效果对比

我们可以对比原始输入图像和基于重建结果的虚拟重光照图像：

def compare_relighting(original_images, albedo, normals, light_directions, index=0):
    """
    对比原始图像和重光照结果
    :param original_images: 原始图像列表
    :param albedo: 反射率图
    :param normals: 法线图
    :param light_directions: 光源方向列表
    :param index: 要对比的图像索引
    """
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
    
    # 显示原始图像
    ax1.imshow(original_images[index], cmap='gray')
    ax1.set_title(f'Original Image {index+1}')
    ax1.axis('off')
    
    # 显示重光照结果
    relit = relight_scene(albedo, normals, light_directions[index])
    ax2.imshow(relit, cmap='gray')
    ax2.set_title(f'Relit Result {index+1}')
    ax2.axis('off')
    
    plt.show()

4. 高级技巧与问题排查

实际应用中会遇到各种问题，以下是常见问题及解决方案：

4.1 处理非Lambertian表面

真实物体往往不完全是漫反射的，可能包含镜面高光。我们可以通过以下方法处理：

1. 鲁棒光度立体法 ：使用更鲁棒的损失函数，减少高光区域的影响
2. 高光检测与去除 ：识别并排除高光像素
3. 多阶段处理 ：先估计漫反射部分，再处理镜面反射

def robust_photometric_stereo(images, light_directions, mask=None, iterations=3):
    """
    鲁棒光度立体法，迭代去除异常值
    """
    k, h, w = images.shape
    if mask is None:
        mask = np.ones((h, w), dtype=bool)
    
    # 初始估计
    normals, albedo = compute_normals_and_albedo(images, light_directions, mask)
    
    for _ in range(iterations):
        # 计算每个图像的重建误差
        errors = []
        for i in range(k):
            relit = relight_scene(albedo, normals, light_directions[i])
            error = np.abs(images[i] - relit)
            errors.append(error)
        errors = np.array(errors)  # [k,h,w]
        
        # 找出误差大的像素（可能是高光）
        median_error = np.median(errors, axis=0)
        outlier_mask = errors > (median_error + 2 * np.std(errors))
        
        # 创建新的图像集，去除异常值
        cleaned_images = images.copy()
        for i in range(k):
            cleaned_images[i][outlier_mask[i]] = np.nan
        
        # 重新计算法线和反射率
        normals, albedo = compute_normals_and_albedo(
            np.nan_to_num(cleaned_images), 
            light_directions, 
            mask
        )
    
    return normals, albedo

4.2 光源方向校准

如果光源方向不够准确，可以使用以下方法改进：

1. 自校准光度立体法 ：同时优化法线和光源方向
2. 使用已知几何的参考物体 ：在场景中放置已知形状的物体来校准光源

def calibrate_light_directions(images, initial_light_directions, mask=None, iterations=5):
    """
    迭代优化光源方向
    """
    k = len(images)
    current_lights = initial_light_directions.copy()
    
    for _ in range(iterations):
        # 固定光源，计算法线和反射率
        normals, albedo = compute_normals_and_albedo(images, current_lights, mask)
        
        # 固定法线，优化光源方向
        for i in range(k):
            # 选择有效像素
            valid = mask & (albedo > 0.1)
            I = images[i][valid]
            N = normals[valid]
            A = albedo[valid]
            
            # 解线性方程组 I = A * (N @ s)
            # 重写为 I = (A * N) @ s
            AN = A[:, np.newaxis] * N
            s = np.linalg.lstsq(AN, I, rcond=None)[0]
            
            # 归一化光源方向
            current_lights[i] = s / (np.linalg.norm(s) + 1e-10)
    
    return current_lights

4.3 从法线重建高度图

有时我们需要从法线图重建实际的3D高度图：

def normals_to_depth(normals, mask=None):
    """
    从法线图重建深度图
    :param normals: [h,w,3] 法线图
    :param mask: [h,w] 可选背景掩码
    :return: [h,w] 深度图
    """
    h, w = normals.shape[:2]
    
    # 计算梯度场
    nz = np.clip(normals[..., 2], 1e-10, 1)  # 避免除以零
    p = -normals[..., 0] / nz  # dz/dx
    q = -normals[..., 1] / nz  # dz/dy
    
    # 初始化深度图
    depth = np.zeros((h, w))
    
    # 从左上角开始积分
    for i in range(1, h):
        depth[i, 0] = depth[i-1, 0] + q[i, 0]
    
    for j in range(1, w):
        depth[0, j] = depth[0, j-1] + p[0, j]
    
    for i in range(1, h):
        for j in range(1, w):
            depth[i, j] = (depth[i-1, j] + q[i, j] + depth[i, j-1] + p[i, j]) / 2
    
    # 归一化深度图
    if mask is not None:
        valid_depths = depth[mask]
        if len(valid_depths) > 0:
            min_d = np.min(valid_depths)
            max_d = np.max(valid_depths)
            depth = (depth - min_d) / (max_d - min_d + 1e-10)
    
    return depth

5. 实际应用案例

让我们看一个完整的应用流程，从原始图像到最终的可视化结果。

5.1 数据准备

假设我们有一个石膏雕像的8张照片，分别从不同方向打光：

# 加载图像
image_files = [f'statue_{i}.jpg' for i in range(1, 9)]
images = [cv2.imread(f, cv2.IMREAD_GRAYSCALE) for f in image_files]
images = np.array(images, dtype=np.float32) / 255.0

# 创建简单背景掩码
avg_intensity = np.mean(images, axis=0)
mask = avg_intensity > 0.1

# 已知光源方向（示例）
light_directions = np.array([
    [0.17, -0.45, 0.88],
    [0.38, -0.38, 0.84],
    [0.55, -0.25, 0.80],
    [0.67, -0.08, 0.74],
    [0.71, 0.12, 0.69],
    [0.64, 0.30, 0.71],
    [0.50, 0.45, 0.74],
    [0.30, 0.55, 0.78]
])
light_directions = light_directions / np.linalg.norm(light_directions, axis=1, keepdims=True)

5.2 计算法线和反射率

# 计算初始法线和反射率
normals, albedo = compute_normals_and_albedo(images, light_directions, mask)

# 鲁棒估计（处理高光）
normals, albedo = robust_photometric_stereo(images, light_directions, mask)

# 后处理
normals = process_normals(normals, albedo, mask)

5.3 结果可视化

# 可视化结果
visualize_normals(normals, mask)
visualize_albedo(albedo, mask)

# 对比原始图像和重光照结果
compare_relighting(images, albedo, normals, light_directions, index=0)
compare_relighting(images, albedo, normals, light_directions, index=4)

# 生成深度图
depth = normals_to_depth(normals, mask)
plt.imshow(depth, cmap='viridis')
plt.title('Depth Map')
plt.colorbar()
plt.show()

5.4 虚拟重光照演示

我们可以交互式地改变光源方向，观察物体外观变化：

from matplotlib.widgets import Slider

fig, ax = plt.subplots()
plt.subplots_adjust(bottom=0.25)

# 初始光源方向
initial_light = np.array([0, 0, 1])
relit = relight_scene(albedo, normals, initial_light)
im = ax.imshow(relit, cmap='gray')

# 创建滑块
ax_slider1 = plt.axes([0.25, 0.1, 0.65, 0.03])
ax_slider2 = plt.axes([0.25, 0.05, 0.65, 0.03])
slider1 = Slider(ax_slider1, 'Light X', -1, 1, valinit=initial_light[0])
slider2 = Slider(ax_slider2, 'Light Y', -1, 1, valinit=initial_light[1])

def update(val):
    x = slider1.val
    y = slider2.val
    z = np.sqrt(1 - x**2 - y**2)  # 保持单位长度
    light = np.array([x, y, z])
    relit = relight_scene(albedo, normals, light)
    im.set_array(relit)
    fig.canvas.draw_idle()

slider1.on_changed(update)
slider2.on_changed(update)

plt.show()