用Python和NumPy手把手实现光度立体法从多张照片重建物体3D表面细节想象一下你只需要用手机环绕物体拍摄几张照片就能自动生成这个物体的3D表面细节——包括每个微小凹凸的法线方向和材质反光特性。这就是光度立体法(Photometric Stereo)的魅力所在。不同于昂贵的3D扫描设备这种方法仅需普通相机和可控光源就能实现惊人的表面细节重建效果。本文将带你用Python和NumPy一步步实现这个计算机视觉中的经典算法。无论你是想为游戏开发快速生成材质贴图还是需要为工业检测分析产品表面缺陷掌握这项技术都能大幅提升你的工作效率。我们将完全从实践角度出发避开复杂的数学推导专注于代码实现和实际问题的解决。1. 准备工作与环境搭建在开始编码前我们需要确保准备好所有必要的工具和数据。光度立体法的核心输入是一组在相同视角、不同光照条件下拍摄的物体照片以及对应的光源方向信息。1.1 安装必要的Python库我们将使用以下Python库来实现算法pip install numpy opencv-python matplotlib scipyNumPy处理所有矩阵运算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, dtypenp.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, axis1, keepdimsTrue)注意光源方向向量应为从物体指向光源的方向且z分量通常为正假设相机沿z轴正方向观察2. 核心算法实现有了图像和光源数据我们现在可以实现光度立体法的核心计算部分。算法主要分为三步计算法线和反射率、后处理、重光照演示。2.1 计算法线和反射率根据Lambertian反射模型图像亮度可以表示为I ρ * (n · s)其中I观察到的像素亮度ρ表面反射率albedon表面法线单位向量s光源方向单位向量对于每个像素我们有多个方程对应多张图像可以解出ρ和n。def compute_normals_and_albedo(images, light_directions, maskNone): 计算法线贴图和反射率图 :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), dtypebool) mask_flat mask.flatten() # [h*w] # 初始化输出 normals np.zeros((3, h*w), dtypenp.float32) albedo np.zeros(h*w, dtypenp.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, axis0) # 计算法线 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, axis2, keepdimsTrue) 1e-10) return normals, albedo2.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], size3) # 2. 只在mask区域内处理 filtered_normals[~mask] normals[~mask] # 3. 重新归一化 norms np.linalg.norm(filtered_normals, axis2) filtered_normals filtered_normals / (norms[..., np.newaxis] 1e-10) return filtered_normals2.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, axis2) shading np.clip(shading, 0, 1) # 限制在[0,1]范围 # 应用反射率 relit albedo * shading return relit3. 结果可视化与分析计算完成后我们需要有效展示和评估结果。法线图通常以RGB颜色表示其中R、G、B通道分别对应法线的x、y、z分量。3.1 可视化法线图import matplotlib.pyplot as plt def visualize_normals(normals, maskNone): 可视化法线图 :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, maskNone): 可视化反射率图 :param albedo: [h,w] 反射率 :param mask: [h,w] 可选背景掩码 plt.imshow(albedo, cmapgray) plt.title(Albedo Map) plt.axis(off) plt.colorbar() plt.show()3.3 重光照效果对比我们可以对比原始输入图像和基于重建结果的虚拟重光照图像def compare_relighting(original_images, albedo, normals, light_directions, index0): 对比原始图像和重光照结果 :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], cmapgray) ax1.set_title(fOriginal Image {index1}) ax1.axis(off) # 显示重光照结果 relit relight_scene(albedo, normals, light_directions[index]) ax2.imshow(relit, cmapgray) ax2.set_title(fRelit Result {index1}) ax2.axis(off) plt.show()4. 高级技巧与问题排查实际应用中会遇到各种问题以下是常见问题及解决方案4.1 处理非Lambertian表面真实物体往往不完全是漫反射的可能包含镜面高光。我们可以通过以下方法处理鲁棒光度立体法使用更鲁棒的损失函数减少高光区域的影响高光检测与去除识别并排除高光像素多阶段处理先估计漫反射部分再处理镜面反射def robust_photometric_stereo(images, light_directions, maskNone, iterations3): 鲁棒光度立体法迭代去除异常值 k, h, w images.shape if mask is None: mask np.ones((h, w), dtypebool) # 初始估计 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, axis0) 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, albedo4.2 光源方向校准如果光源方向不够准确可以使用以下方法改进自校准光度立体法同时优化法线和光源方向使用已知几何的参考物体在场景中放置已知形状的物体来校准光源def calibrate_light_directions(images, initial_light_directions, maskNone, iterations5): 迭代优化光源方向 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, rcondNone)[0] # 归一化光源方向 current_lights[i] s / (np.linalg.norm(s) 1e-10) return current_lights4.3 从法线重建高度图有时我们需要从法线图重建实际的3D高度图def normals_to_depth(normals, maskNone): 从法线图重建深度图 :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 depth5. 实际应用案例让我们看一个完整的应用流程从原始图像到最终的可视化结果。5.1 数据准备假设我们有一个石膏雕像的8张照片分别从不同方向打光# 加载图像 image_files [fstatue_{i}.jpg for i in range(1, 9)] images [cv2.imread(f, cv2.IMREAD_GRAYSCALE) for f in image_files] images np.array(images, dtypenp.float32) / 255.0 # 创建简单背景掩码 avg_intensity np.mean(images, axis0) 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, axis1, keepdimsTrue)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, index0) compare_relighting(images, albedo, normals, light_directions, index4) # 生成深度图 depth normals_to_depth(normals, mask) plt.imshow(depth, cmapviridis) plt.title(Depth Map) plt.colorbar() plt.show()5.4 虚拟重光照演示我们可以交互式地改变光源方向观察物体外观变化from matplotlib.widgets import Slider fig, ax plt.subplots() plt.subplots_adjust(bottom0.25) # 初始光源方向 initial_light np.array([0, 0, 1]) relit relight_scene(albedo, normals, initial_light) im ax.imshow(relit, cmapgray) # 创建滑块 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, valinitinitial_light[0]) slider2 Slider(ax_slider2, Light Y, -1, 1, valinitinitial_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()
