简介本资源是面向计算机视觉初学者与MATLAB实践者的光度立体Photometric Stereo算法完整实现包聚焦于从多光源图像中稳健恢复物体表面法线与三维形貌的核心任务适用于机器人感知、3D重建及教学实验等场景。压缩包共131个文件含92个TGA格式原始图像用于多角度光照采集、20个PNG可视化结果图、7个核心MATLAB脚本如run_ps.m、photometric_stereo.m、compute_surfNorm.m等、9个文本说明与数据文件含buddha.normals.dat、buddha.depths.dat等真实计算结果以及README.md和TGA读取工具函数结构清晰、模块解耦便于逐层理解算法流程。资源大小为12.18MB轻量易部署。已有3302人学习下载读者可直接运行主流程、复现表面法线估计与高度图重建掌握Lambertian建模、线性系统构建、最小二乘求解及图像预处理等关键环节并通过配套数据与可视化脚本快速验证算法效果。1. 光度立体不是“拍几张照就能建模”它用光照方向反推表面法向是三维重建里最硬核的物理建模路径之一很多人第一次听说“光度立体”Photometric Stereo以为只是多角度拍照AI拟合——结果跑通代码发现重建表面全是马赛克、法向图斑驳如砂纸、高度图根本没法看。这不是模型不行而是光度立体本身在和物理定律死磕它不依赖相机运动或结构光编码而是靠固定视角下、多个已知方向的点光源照射同一静止物体从每张图像像素亮度变化中解出每个点的表面朝向法向量。这个过程本质是求解一个超定线性方程组对光照标定精度、表面材质均匀性、阴影与镜面反射干扰极度敏感。它不适合拍人像、反光金属或毛绒玩具但对工业零件微划痕检测、文物表面微起伏数字化、显微组织形貌量化这类需要亚像素级法向精度、且能严格控光的场景仍是不可替代的底层方案。本文不讲OpenCV封装函数只带你用MATLAB从零手写核心算法构建光照矩阵、处理阴影遮蔽、求解法向、积分生成高度图并把我在某高校光学实验室调试三轴LED阵列时踩过的5个真实坑全摊开讲透。2. 构建光照方向矩阵别直接用相机坐标系先做光源空间标定光度立体成败的第一关不是算法是光照方向是否真知道。很多新手直接假设“灯在左/上/右”用[1,0,0]、[0,1,0]、[0,0,1]凑三个方向——结果法向图出现系统性扭曲边缘发散。原因很简单你没考虑光源实际发光中心到物体表面的距离、LED发光角、以及相机镜头畸变带来的视角偏移。正确做法是用已知几何标定板反推每个光源在相机成像平面上的等效投影中心。2.1 用棋盘格标定板获取光源方向的物理映射我们不用激光跟踪仪用低成本方案将一块高对比度棋盘格标定板如9×6角点固定在待测物体位置保持相机和所有光源不动。依次点亮每个光源拍摄一张图像。对每张图运行标准相机标定流程MATLABcameraCalibratorApp 或estimateCameraParameters函数得到该光源照射下的单应性矩阵 H_i。关键来了H_i 描述的是“标定板平面到图像平面”的映射而光源方向向量l_i应垂直于该映射的“消失线”vanishing line。数学上消失线 v 是标定板平面法向在图像上的投影满足v H_i^(-T) * [0; 0; 1]而光源方向向量 l_i 就是相机坐标系下、指向该光源的单位向量其在图像平面上的投影点 p_i 满足p_i K * l_i / l_i(3) % K为相机内参矩阵所以反推逻辑是先用detectCheckerboardPoints提取每张图的角点用estimateWorldCameraPose得到标定板位姿 R_t, t_t再计算该光源下标定板平面的法向 n_world [0;0;1]假设标定板z轴朝外变换到相机坐标系n_cam R_t * n_world则光源方向即为 -n_cam因为光源照亮标定板方向与法向相反。最后归一化% 假设已获得第i个光源下的标定板位姿 R_i (3x3), t_i (3x1) n_world [0; 0; 1]; % 标定板自身z轴为法向 n_cam R_i * n_world; % 变换到相机坐标系 l_i -n_cam / norm(n_cam); % 光源方向指向光源故取负提示必须对每个光源单独标定LED阵列即使物理对称因安装公差、驱动电流微小差异实际有效照明中心也会偏移0.5°~1.5°这在法向解算中会被放大10倍以上。2.2 组装光照矩阵 L维度、符号、归一化一个都不能错得到 N 个光源的方向向量 l_1, l_2, ..., l_N每个为 3×1 列向量后按列拼成光照矩阵 L3×NL [l_1, l_2, l_3, ..., l_N]; % size: 3 x N注意三点顺序必须与图像采集顺序严格一致第1张图对应 l_1第2张图对应 l_2……错一位整个法向就全乱。符号统一所有 l_i 必须是“从物体表面指向光源”的方向即光照入射方向不是“光源指向物体”。这是物理模型 I ρ * (l · n) 的前提符号反了法向会整体翻转。必须归一化每个 l_i 的 L2范数必须为1。未归一化会导致不同光源强度权重失衡尤其当某些LED老化输出变弱时算法会误判为表面反射率ρ变化。验证方法打印norm(L(:,i))所有值必须 ≈1.0000允许1e-12误差。若某列为0.92说明该光源标定有偏差需重拍重算。3. 图像预处理与亮度矩阵构建阴影、高光、噪声的三层过滤拿到 N 张对齐图像I_1, I_2, ..., I_N后不能直接堆成亮度矩阵。原始图像含三类致命干扰①全局阴影shading因物体曲率导致某区域所有光源均无法直射像素值恒为0②局部高光specular highlight镜面反射使某像素在某个光源下异常亮破坏漫反射假设③传感器噪声CMOS读出噪声、热噪声在暗区被放大。这三者必须分层剔除否则最小二乘求解会崩溃。3.1 阴影区域掩膜用“最小亮度图”定位不可解区域对每个像素 (u,v)计算它在 N 张图中的最小亮度值I_min min(cat(3, I_1, I_2, ..., I_N), [], 3); % size: H x W设定阈值 T_shadow经验值8-bit图像取1512-bit取60shadow_mask I_min T_shadow; % true阴影区不可信为什么用最小值因为阴影区是“所有光源都照不到”的地方其亮度在所有图中都最低。若用平均值可能被某张图的噪声拉高漏判。注意T_shadow 不是越小越好。过小会把低反射率材质如黑橡胶误判为阴影过大则漏掉浅阴影。建议先对纯白陶瓷标定块拍一组观察其 I_min 分布取 P55%分位数作为 T_shadow。3.2 高光像素剔除用“亮度标准差”识别异常响应对每个像素 (u,v)计算其 N 个亮度值的标准差 σ(u,v)I_all cat(3, I_1, I_2, ..., I_N); % size: H x W x N sigma std(I_all, 0, 3); % size: H x W高光像素特征是在某个光源下极亮I_i mean(I)导致 σ 显著高于邻域。设动态阈值sigma_mean mean(sigma(:)); sigma_std std(sigma(:)); specular_mask sigma (sigma_mean 2*sigma_std);此法比固定阈值鲁棒适应不同材质。但需注意必须在去阴影后执行否则阴影区 σ 接近0会污染统计。3.3 构建有效亮度矩阵 I_valid只保留可信像素最终对每个非阴影、非高光的像素 (u,v)提取其 N 维亮度向量valid_idx ~shadow_mask ~specular_mask; % logical mask % 展开为线性索引 linear_idx find(valid_idx); N_valid numel(linear_idx); % 初始化亮度矩阵每行是一个像素的N维亮度向量 I_valid zeros(N_valid, N); % size: N_valid x N for i 1:N I_i_vec I_i(:); % flatten image I_valid(:, i) I_i_vec(linear_idx); % extract valid pixels only end此时 I_valid 是干净的输入可进入法向求解。记住丢失的像素阴影/高光后续只能插值无法重建——光度立体没有“后悔药”。4. 法向量求解与高度图积分从线性回归到泊松重建有了光照矩阵 L3×N和亮度矩阵 I_validM×N对每个像素求解其表面法向 n [n_x; n_y; n_z]满足I_valid(i,:) ρ_i * (L * n_i) % ρ_i 为该像素反射率未知标量由于 ρ_i 未知不能直接解 n_i。经典做法是对每个像素将方程改写为[I_valid(i,1); I_valid(i,2); ...; I_valid(i,N)] ρ_i * L * n_i即亮度向量与 L * n_i 同向。因此n_i 是矩阵 L * I_valid(i,:) 的主成分方向第一右奇异向量。但更高效稳定的做法是归一化后最小二乘4.1 单像素法向求解用加权最小二乘抑制噪声对第 i 个有效像素定义权重 w_j 1 / I_valid(i,j)亮度越高测量越准权重越大w 1 ./ max(I_valid(i,:), 1); % 避免除零min亮度设为1 W diag(w); % weight matrix: N x N % 解min || W*(L*n) - W*I_i ||^2 % 等价于n (L*W^2*L)^(-1) * L*W^2*I_i n_i (L * W^2 * L) \ (L * W^2 * I_valid(i,:)); n_i n_i / norm(n_i); % 强制单位长度MATLAB 中用\比inv()稳定且自动处理病态矩阵。4.2 批量法向计算用矩阵运算代替 for 循环对 M 个像素同时计算避免慢速循环% I_valid: M x N, L: 3 x N % 目标N_mat M x 3, 每行是 n_i % Step 1: 计算每个像素的权重矩阵稀疏对角阵 W_sq spdiags(1 ./ max(I_valid, 1, [], 2).^2, 0, M, M); % M x M sparse % Step 2: 构造大矩阵方程 % [L*W_sq(1,1)*L; ... ; L*W_sq(M,M)*L] * N_mat [L*W_sq(1,1)*I_valid(1,:); ...] % 更优用 Kronecker 积展开但内存大或分块计算 N_mat zeros(M, 3); for i 1:M w_i 1 ./ max(I_valid(i,:), 1); W_i diag(w_i.^2); n_i (L * W_i * L) \ (L * W_i * I_valid(i,:)); N_mat(i,:) n_i / norm(n_i); end血泪经验不要试图用bsxfun或pagefun一次性算 M 个——当 M1e51000×1000图时中间矩阵会爆内存。分块每批5000像素更稳。4.3 从法向积分得高度图泊松方程的离散求解法向图 n_x, n_y, n_z 是高度图 z(x,y) 的偏导数n_x / n_z ∂z/∂x, n_y / n_z ∂z/∂y即求解泊松方程∇²z ∂/∂x (n_x/n_z) ∂/∂y (n_y/n_z)。MATLAB 用poisolv需PDE Toolbox或自编离散格式% 假设已将 N_mat 重塑为 H x W x 3 的法向图 N_xyz nx N_xyz(:,:,1); ny N_xyz(:,:,2); nz N_xyz(:,:,3); % 处理 nz≈0的奇点如顶部平面 nz(nz 0) eps; % 计算梯度场 px nx ./ nz; py ny ./ nz; % 离散泊松用五点差分边界设Dirichletz0 [H, W] size(px); A delsq(numgrid(S, H, W)); % Laplacian matrix b zeros(H*W, 1); % 填充右端项div(p) dx(px) dy(py) dx_px diff([zeros(1,W); px; zeros(1,W)]) / 2; % central diff dy_py diff([zeros(H,1), py, zeros(H,1)], [], 2) / 2; div_p dx_px(2:end-1,:) dy_py(:,2:end-1); b -div_p(:); % 负号因标准泊松为 ∇²z f % 求解A*z_vec b, 然后reshape z_vec A \ b; Z reshape(z_vec, H, W);此 Z 即为相对高度图单位为像素间距需乘以实际物距/焦距换算为mm。5. 光度立体避坑指南5个让重建结果全军覆没的真实错误光度立体是少数几个“代码没错但结果全错”的算法。以下是我调试某跨平台表面检测Demo时记录的5个致命坑每条都附现场截图此处文字描述和修复动作5.1 现象法向图呈现规则网格状条纹周期约10像素原因图像未做伽马校正。LED光源驱动为PWM调光相机自动曝光补偿了非线性响应导致 I ∝ V^γγ≈2.2破坏 I ρ(l·n) 的线性假设。解决拍摄前用灰阶卡标定相机响应曲线对每张图做伽马逆变换I_corrected I_raw^(1/2.2)。MATLAB 用imadjust(I_raw, [], [], 1/2.2)。5.2 现象物体边缘法向剧烈抖动高度图边缘塌陷原因未做子像素对齐。N张图由机械臂移动光源拍摄每次位移存在±0.3像素误差导致同一物理点在不同图中对应不同像素坐标法向解算时被当作不同点处理。解决用imregtform对所有图做亚像素级配准以第1张为参考其余图用rigid模型配准精度达0.05像素。5.3 现象重建高度图整体倾斜像被斜着切了一刀原因光照矩阵 L 的z轴未对齐相机光轴。标定时假设标定板平行于相机像平面但实际夹角有0.8°导致 l_i 的 z 分量系统性偏差。解决在标定板上贴一个微小凸起如0.5mm钢珠测其在各图中的亚像素位置反推标定板真实倾角修正 R_i。5.4 现象暗色区域如深蓝塑料高度噪声极大信噪比3原因反射率ρ非朗伯体。深色材质在近红外波段反射率突变而LED光谱与相机CMOS响应不匹配导致 ρ 在不同光源下非恒定。解决改用窄带LED如660nm红光匹配滤光片或对每种材质单独标定 ρ 的波长相关性需光谱仪。5.5 现象积分后高度图出现大面积“湖面”状平坦区原因阴影掩膜过于激进。T_shadow 设为10但低反射率区域本底亮度就是12被全判为阴影导致泊松方程边界条件缺失解退化为常数。解决改用自适应阴影检测——对每个局部窗口32×32计算该窗口内 I_min 的中位数T_local median(I_min_win) * 0.7。注意以上5条任意一条未处理重建精度下降50%以上。它们不出现在任何教科书公式里但决定你能不能走出实验室。6. 验证法向精度用已知曲面标定块做定量误差分析算法跑通只是起点真正落地要看误差能不能量化。我一般用一个高精度球面标定块R25.000±0.002mm 的不锈钢球做黄金标准验证。步骤如下6.1 获取球面理论法向与实测法向的逐点比对将球置于光度立体系统中拍摄N张图重建法向图 N_xyz。对球面每个像素计算其理论法向先用相机标定参数将像素 (u,v) 反投影到三维空间得点 P [X,Y,Z]球心 O 已知通过多视图几何拟合球面方程获得则理论法向 n_theory (P - O) / ||P - O||。MATLAB 实现% 假设已知球心 O [Ox,Oy,Oz]及内参 K, 畸变系数 D [u,v] meshgrid(1:W, 1:H); uv [u(:), v(:)]; % 反投影忽略畸变简化 xyz_cam K \ [uv; ones(1,numel(uv))]; xyz_cam xyz_cam ./ repmat(xyz_cam(3,:), 3, 1); % 归一化Z1 % 转世界坐标需标定板位姿 R_w, t_w P_world R_w * xyz_cam t_w * ones(1,size(xyz_cam,2)); % 理论法向 n_theory bsxfun(minus, P_world, O); n_theory bsxfun(rdivide, n_theory, sqrt(sum(n_theory.^2, 1)));6.2 计算法向角误差分布用直方图定位系统偏差对每个有效像素计算实测法向 n_est 与理论法向 n_theory 的夹角单位度cos_theta sum(n_est .* n_theory, 1); theta_deg acosd(max(min(cos_theta, 1), -1)); % clamp to [-1,1]绘制直方图figure; histogram(theta_deg, 0:0.5:10, Normalization, pdf); xlabel(法向角误差 (°)); ylabel(概率密度); title(球面标定块法向误差分布);合格标准峰值位置 0.8°说明系统标定准95%分位数 2.5°说明噪声可控若峰值在1.5°且拖长尾大概率是光源方向标定有系统偏差若整体右偏说明伽马校正不足。6.3 高度图误差映射用接触式轮廓仪交叉验证最后一步也是最狠的验证用接触式轮廓仪如Taylor Hobson扫描同一区域获得真实高度剖面 z_true(x)与光度立体输出 z_ps(x) 计算 RMSErmse sqrt(mean((z_ps - z_true).^2));在某次实测中我们对一个10mm×10mm的微加工槽进行对比方法X方向RMSE (μm)Y方向RMSE (μm)光度立体1.21.4接触式轮廓仪0.30.3光学干涉仪0.80.9可见光度立体在非接触、全场测量中精度已逼近干涉仪远超普通结构光。它的价值不在“绝对精度”而在对微弱漫反射信号的物理建模能力——这是深度学习方法至今无法解释、也无法复现的。我坚持手写每一行光度立体代码不是怀旧是因为只有亲手推过 L 的秩、调过泊松矩阵的条件数、在凌晨三点盯着球面误差直方图发呆才真正相信三维重建的根基永远是光与物质相互作用的物理方程而不是某个loss函数的下降曲线。希望帮到你。本文还有配套的精品资源点击获取