当前位置：首页 > news >正文

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

news 2026/5/24 14:57:12

用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 光源方向测量

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

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

假设我们已经测得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表面

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

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

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 光源方向校准

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

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

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

查看全文

http://www.jsqmd.com/news/877920/