简介：这是一份基于C++与OpenCV实现的光度立体算法可运行源码，面向计算机视觉初学者和三维重建研究者，用于通过多角度光照图像恢复物体表面法向量、计算广义反射率与梯度图，解决从算法原理到工程实现的落地问题。压缩包共12个文件、约69KB，包含核心C++源文件、Makefile构建脚本、7张测试输入与结果图像（法向量图、梯度图、反射率图）及项目配置文件，整体结构精简，便于快速编译和运行验证。源码覆盖球系坐标转笛卡尔坐标、光照方向矩阵初始化与归一化、线性方程组求解等关键环节，并配有完整的可执行示例和对应输出图，读者可直接观察每个步骤的效果。已有165人学习浏览，适合需要从数学推导过渡到OpenCV编码实践的开发者，既可在实际工程中参考调试，也可作为课程设计与课题研究的基础代码。
1. 光度立体算法实现：一套 OpenCV C++ 可运行源码，把多角度光照拍出三维法线
光度立体这个词听着唬人，拆开其实就一件事：让相机不动，换 N 个方向的光源拍 N 张照片，从这几张图里解出每个像素的表面法线方向和反照率。相比双目、结构光，它不需要第二台相机，也不需要投影设备，所以很多做表面缺陷检测、文物理貌复原、材质采集的从业者都会先试这条路。这里整理的光度立体算法实现，是一套基于 OpenCV C++ 的可运行源码，输入一组多光源图像和光源方向文件，直接输出法线贴图和反照率图。适合刚接触光度立体的同学照着跑通流程，也适合需要把算法嵌进现有视觉检测流水线的工程师直接改参数。
2. 光度立体原理：从朗伯体模型到最小二乘求解，关键在光源方向矩阵
2.1 朗伯体假设与图像辐照度方程
光度立体能成立，前提是目标表面近似朗伯体，也就是表面把入射光均匀散射到各个方向，观察方向不影响亮度。这样每个像素的灰度只取决于三个量：表面反照率 ρ、光源强度 E、以及光源方向 l 与表面法线 n 的夹角余弦。写成一个基础方程就是：
I = ρ · E · max(l·n, 0)
这里的 l 是从物体表面指向光源的单位向量，n 是表面的单位法线。max 的意义是光源照不到背面时，亮度直接为 0。对大多数陶瓷、塑料、未抛光的金属件来说，这个假设在入射角不太大的时候是够用的。
单张照片没法把 ρ 和 n 拆开，因为一个像素只有一个亮度值，却有两个未知数。光度立体把光源方向换 N 次，等于给每个像素测出了 N 个亮度值。只要 N≥3 且光源方向不共面，就能组成一个超定方程组，把 ρ 和 n 解出来。实际工程里我一般取 8 到 12 张，光方向尽量均匀分布在以物体为中心的上半球，避免所有光源挤在同一个平面附近。
2.2 多光源采集与矩阵求解形式
如果每个像素单独看，可以写成线性方程 I = L·m。这里的 L 是 N 行 3 列的矩阵，每一行是一个光源方向；I 是 N 行 1 列的向量，是同一个像素在 N 张图里的灰度；未知量 m 是一个三维向量，它等于 ρ·n。因为 n 是单位向量，所以求解 m 之后，法线就是 m 除以它的模长，反照率就是 m 的模长。
L 不是方阵，不能直接求逆，所以要解这个超定方程组。常见做法是最小二乘，对应正规方程为 m = (L^T L)^(−1) L^T I。正规方程算起来简单，但 L 的条件数不好时，数值上容易出现小奇异值被放大。源码里用 SVD 分解解同样的最小二乘问题，虽然速度比正规方程慢一点，但稳定性好很多，处理 8 到 12 个光源的这种小矩阵完全不是瓶颈。
2.3 为什么这套源码选 OpenCV 和 C++
同类开源实现很多是用 Python 写的，但实际落到工业场景时，C++ 路线的价值主要体现在三件事。第一，OpenCV 的 cv::solve 支持多列同时求解，一次调用就把整张图像的线性系统算完；第二，图像读取、灰度化、掩膜、保存 PNG 这些 IO 操作都是现成的，不用额外引入其他依赖；第三，C++ 写的算法可以直接封装成动态库，挂到现有视觉软件里，不用在一套系统里再塞一个 Python 运行时。
这套源码没有依赖 OpenCV 之外的三方库，编译时只要保证 OpenCV 3.x 或 4.x 可用。CPU 上处理 1280×1024 的图像，8 张输入，Release 模式下整条流程跑完大概一两秒，主要耗在内存拷贝和法线图写入上，计算本身很快。
3. 可运行源码结构：数据组织、模块划分与编译方式
3.1 输入数据长什么样：多光源图像集与光源方向文件
源码的输入是三类文件：一组同一视角、不同光源方向拍摄的灰度图或彩色图，一个光源方向文本文件，以及一张可选的掩膜图像。图像用普通 PNG 或 JPG 都行，但强烈建议用 PNG 或直接保存的 TIFF，JPG 压缩在高光和暗部容易引入噪点，影响法线精度。
光源方向文件的格式很简单，每行三个浮点数，分别代表 x、y、z，一行对应一张图像：
0.00 0.00 1.00 0.35 0.00 0.94 -0.35 0.00 0.94 0.00 0.35 0.94 0.00 -0.35 0.94这些方向都必须是从物体指向光源的单位向量。源码在读入时会做一次归一化，所以文件里写没写单位长度都无所谓，但方向正负绝对不能反。
掩膜图像的作用是排除背景、高光或阴影区域。掩膜里像素值为 0 的位置不参与求解，法线图对应位置填默认的 (0,0,1)。如果你没有掩膜，可以传一个全 255 的白色图片，或者让源码通过暗部阈值自动生成。
3.2 源码目录与模块划分
源码包解压后结构大概是这样：
photometric_stereo/ ├─ data/ │ ├─ sample/ │ │ ├─ light_00.png │ │ ├─ light_01.png │ │ └─ lights.txt │ └─ sphere/ │ ├─ light_00.png │ └─ lights.txt ├─ include/ │ └─ photometric_stereo.h ├─ src/ │ ├─ main.cpp │ ├─ load_data.cpp │ └─ photometric_stereo.cpp ├─ scripts/ │ └─ make_sphere_dataset.py ├─ CMakeLists.txt └─ README.mdmain.cpp 负责参数解析和整体流程调度；load_data.cpp 负责读取图像、光源方向、掩膜；photometric_stereo.cpp 是核心算法，包含矩阵构建、SVD 求解、法线和反照率计算。scripts 里放了合成球体数据集脚本，用来验证算法正确性，这个脚本在文章最后一节会重点讲。
3.3 编译方式与 CMake 配置
源码用 CMake 构建。OpenCV 3.x 和 4.x 都兼容，构建时需要指定 OpenCV 的安装路径，或者让 CMake 通过系统环境变量自动找到。CMakeLists.txt 的核心内容如下：
cmake_minimum_required(VERSION 3.10) project(photometric_stereo) set(CMAKE_CXX_STANDARD 11) set(CMAKE_CXX_STANDARD_REQUIRED ON) find_package(OpenCV REQUIRED) add_executable(photo_stereo src/main.cpp src/load_data.cpp src/photometric_stereo.cpp ) target_include_directories(photo_stereo PRIVATE include) target_link_libraries(photo_stereo PRIVATE ${OpenCV_LIBS})find_package 会去查找 OpenCVConfig.cmake。如果你的 OpenCV 安装位置特殊，编译前需要设置 OPENCV_DIR 环境变量指向 OpenCV 的 cmake 目录。设置好后，Linux 或 macOS 上直接编译：
mkdir build && cd build cmake -DCMAKE_BUILD_TYPE=Release .. cmake --build . -j$(nproc) ./photo_stereo "../data/sample/*.png" ../data/sample/lights.txt ../data/sample/mask.png ../results命令行第一个参数是图像路径的通配符，必须用引号包住；第二个是光源方向文件；第三个是掩膜路径；第四个是输出目录。程序内部会用 cv::glob 展开通配符，并对文件名排序，保证图像顺序和光源文件行顺序一一对应。如果你发现法线图是乱的，大概率就是这里顺序错位了。
4. 核心实现拆解：批量 SVD 求解与法线贴图输出
4.1 图像读取、灰度化与掩膜预处理
主流程第一步是把所有输入图像读进来，统一转为灰度图，并转换到 CV_64F 浮点，值域压到 0 到 1 之间。这样做是为了让后续求解的数值范围稳定。代码如下：
std::vector<cv::Mat> readImages(const std::vector<std::string>& paths) { std::vector<cv::Mat> images; for (const auto& p : paths) { cv::Mat img = cv::imread(p, cv::IMREAD_GRAYSCALE); if (img.empty()) { throw std::runtime_error("failed to load image: " + p); } img.convertTo(img, CV_64F, 1.0 / 255.0); images.push_back(img); } return images; }这里用 IMREAD_GRAYSCALE 直接读成单通道，省掉一次 cvtColor 调用。如果你的输入是 16 位 TIFF，需要把缩放系数从 1.0/255.0 改成 1.0/65535.0。很多初次跑的人会在这里翻车：图像看起来正常，但亮度值整体偏大或偏小，导致法线偏向某一个方向。读成单通道后，整个求解都基于灰度亮度，不用关心彩色通道。
掩膜读取更简单，灰度图转成 CV_8U，然后调用 cv::compare 把非零区域变成掩膜。如果你没有掩膜，就给一张全白图，这样所有像素都参与求解。注意掩膜必须和输入图像尺寸一致，否则后续矩阵列数对不上，程序会在运行时直接抛异常。
4.2 构建光度矩阵并批量求解法线
这是整个源码最核心的部分。先把 N 张图拉成一维行向量，堆成 N 行 M 列的矩阵 I，其中 M 是像素总数。光源方向堆成 N 行 3 列的矩阵 L。然后调用 cv::solve 一次求出所有像素的 G，G 是 3 行 M 列，每一列就是该像素的 ρ·n。
void solveNormalMap(const std::vector<cv::Mat>& images, const std::vector<cv::Vec3d>& lights, const cv::Mat& mask, cv::Mat& normalMap, cv::Mat& albedoMap) { int rows = images[0].rows; int cols = images[0].cols; int pixelCount = rows * cols; int lightCount = (int)images.size(); cv::Mat I(lightCount, pixelCount, CV_64F); for (int i = 0; i < lightCount; ++i) { cv::Mat row = I.row(i); images[i].reshape(1, 1).copyTo(row); } cv::Mat L(lightCount, 3, CV_64F); for (int i = 0; i < lightCount; ++i) { L.at<double>(i, 0) = lights[i][0]; L.at<double>(i, 1) = lights[i][1]; L.at<double>(i, 2) = lights[i][2]; } cv::Mat G; cv::solve(L, I, G, cv::DECOMP_SVD); normalMap.create(rows, cols, CV_32FC3); albedoMap.create(rows, cols, CV_64F); for (int idx = 0; idx < pixelCount; ++idx) { int r = idx / cols; int c = idx % cols; if (mask.at<uchar>(r, c) == 0) { normalMap.at<cv::Vec3f>(r, c) = cv::Vec3f(0.0f, 0.0f, 1.0f); albedoMap.at<double>(r, c) = 0.0; continue; } double gx = G.at<double>(0, idx); double gy = G.at<double>(1, idx); double gz = G.at<double>(2, idx); double len = std::sqrt(gx * gx + gy * gy + gz * gz); albedoMap.at<double>(r, c) = len; if (len > 1e-12) { normalMap.at<cv::Vec3f>(r, c) = cv::Vec3f( (float)(gx / len), (float)(gy / len), (float)(gz / len)); } else { normalMap.at<cv::Vec3f>(r, c) = cv::Vec3f(0.0f, 0.0f, 1.0f); } } }这段有三个要点。第一，images[i].reshape(1, 1).copyTo(row) 的作用是把二维灰度图展平成一维，直接放进 I 的第 i 行，省掉手写双层 for 循环。第二，cv::solve 的 DECOMP_SVD 表示用 SVD 求解最小二乘，当 L 接近奇异时，SVD 仍能给出一个稳定的数值解，不会因为正规方程里的矩阵求逆而崩掉。第三，len 小于 1e-12 的像素认为是无效像素，填默认法线，避免除零产生 NaN。
在性能方面，真正耗时的不是 cv::solve 本身，而是后面那个逐像素 for 循环。如果分辨率太高，可以考虑把图像分成 512×512 的块，每块单独建 I 矩阵、单独求解，内存占用会低很多。源码默认一次性处理全图，4K 分辨率、8 张图的内存占用大约在 1GB 左右，这个量级大多数开发机都扛得住。
4.3 法线结果转 RGB 法线贴图与反照率图
法线是三维向量，每个分量范围在 -1 到 1 之间，直接保存成普通 8 位 PNG 会把负值截掉。需要先把分量映射到 0 到 255。代码如下：
cv::Mat visualizeNormalMap(const cv::Mat& normalMap) { cv::Mat vis(normalMap.size(), CV_8UC3); for (int r = 0; r < normalMap.rows; ++r) { for (int c = 0; c < normalMap.cols; ++c) { cv::Vec3f n = normalMap.at<cv::Vec3f>(r, c); cv::Vec3b out; out[0] = cv::saturate_cast<uchar>((n[2] + 1.0) * 127.5); out[1] = cv::saturate_cast<uchar>((n[1] + 1.0) * 127.5); out[2] = cv::saturate_cast<uchar>((n[0] + 1.0) * 127.5); vis.at<cv::Vec3b>(r, c) = out; } } return vis; }OpenCV 里 Mat 的通道顺序是 BGR，所以这里 out[0] 对应蓝色通道，保存的是法线 z 分量。这样做的结果是：正对相机的平面法线是 (0,0,1)，映射后蓝色通道是 255，看起来就是一片蓝色，这是法线贴图的标准观感。out[2] 对应红色通道，保存的是 x 分量。
如果你要用 Blender、Unity 这类工具读取法线贴图，它们内部默认的通道约定可能是 R=x、G=y、B=z，那就要做一次通道重排。源码里把原始浮点法线也保存了一份，避免 PNG 可视化过程中的精度损失。保存使用 OpenCV 的 EXR 格式：
cv::imwrite("normal_map.exr", normalMap);前提是你的 OpenCV 编译时带有 OpenEXR 支持。大多数通过包管理器安装的版本都带。如果不带，可以用 cv::FileStorage 写成 YAML，或者直接写一个自定义的二进制文件，把 float 数据原样落盘。
5. 避坑指南：光度立体常见问题与排查清单
5.1 现象：法线方向整体颠倒，反照率图变“负片”
如果你跑出来的法线贴图整体颜色反了，球面边缘的凸起变成凹陷，先别怀疑算法，八成是光源方向的正负约定写反了。光度立体的约定是从物体表面指向光源。如果你的光源文件里存的向量方向是从光源指向物体，那么求解出来的法线整体会反向。
解决办法是检查 lights.txt 的第一行。假设第一张图是从正前方打光，那么光源方向应该是 (0,0,1)，而不是 (0,0,-1)。如果方向文件格式没错，那就要看读入代码里有没有多做一次翻转。最简单的方法是拿一个已知形状的粗糙球体跑一遍，球体的凸面必须外凸，如果变成内凹，就把所有光源方向乘个 -1。
5.2 现象：法线贴图出现大量颗粒状噪点
这种问题在高光物体上最常见。表面越接近镜面，离开朗伯体假设就越远。当光源方向恰好让高光进入相机时，那个像素的亮度远大于朗伯体预测值，代入线性系统后会把法线带偏。另外，相机如果用 sRGB 或者开了自动增益，图像亮度与光照能量根本不是线性关系，也会导致同样的噪点。
想解决，第一步是确认输入图像有没有经过 gamma 或自动白平衡处理。理想输入是线性灰度图，最好是原始 RAW 转 TIFF，不要用 JPG。第二步是生成掩膜，把高光区域排除掉。可以先用一个固定阈值把亮度大于 250 的像素标记为无效。第三步是改用多次曝光，保证暗部不欠曝、亮部不过曝。
我在实际处理陶瓷表面时有一个经验：不要只靠阈值去高光，因为朗伯体亮部和真正高光的灰度值可能有重叠。更好的是多拍一组交叉偏振光，把镜面反射分量物理上滤掉。但这套基础源码不做偏振输入，所以工程上先用掩膜规避。
5.3 现象：表面凹凸像被“压扁”或局部反向
如果法线方向从整体上看是合理的，但某些区域凹凸反转，通常是光源方向标定不精确。光源方向每偏几度，恢复出的法线就会跟着偏，尤其是在入射角比较大的边缘区域。使用手搭的光源架子，最容易出现的问题是光源实际位置和理论位置不符。
解决思路有两种。第一种是用一个镜面球做标定，从图像里的高光位置反推每个光源的方向。第二种是减少入射角过大的光源，控制在 60 度以内。入视角越大，法线误差对光源方向误差越敏感。还有一个容易被忽略的点：所有光源方向必须定义在同一个坐标系里。如果光源方向是在相机坐标系下标的，而掩膜或法线图用图像坐标系，就会混出一个畸形的结果。源码默认光源方向在相机坐标系下，输出法线也在这个坐标系。
5.4 现象：掩膜把所有像素都排除了
掩膜读取后，参与求解的像素数量为 0，程序没有报错，但输出是一张全蓝的默认法线图。这种问题出现在自动生成掩膜的场景。如果用“图像灰度最小值 < 阈值”来判断有效像素，阈值设太高，比如 50，很多正常阴影区域会被当成无效像素。
解决方法是先把掩膜可视化出来，确认是不是你想保留的区域。调低阈值到 5 或 10，再跑一次。另外注意，掩膜图像的值域是 0 到 255，读取后不要直接拿像素值做逻辑判断，要用 255 代表有效。有的版本习惯写成 255 有效，有的写成 0 有效，源码里固定为 0 是排除、非 0 是保留，这个在 README 里写得很清楚。
5.5 现象：Debug 模式运行极慢，Release 模式内存偏高
这个问题不算算法问题，但很影响体验。Debug 模式下 cv::solve 和逐像素 for 循环都会慢好几倍，1280×1024 的图像可能要等十几秒。Release 模式很快，但内存占用会突然涨上去。原因在于读入的图像全部转成了 CV_64F，一张 4K 灰度图的内存是 8.3M×8 字节，约 64MB，8 张就是 512MB，再算上 I 矩阵和法线图，整体占用轻松超过 1GB。
解决办法有两个。一是用 Release 编译运行，这应该是最基本的要求。二是如果内存仍然吃紧，把图像降采样到 1024 宽再跑，法线图虽然细节少一点，但整体趋势完全够用。对算法验证来说，降采样是更优先的选择。
6. 一个验证技巧：先用合成球体跑通全套管线，再碰真实数据
6.1 为什么先用合成球体
真实场景里的每个误差源都会叠加：光源方向标定误差、相机非线性、高光、阴影、环境光。如果一上来就用真实物体调算法，出了问题根本分不清是哪一环。合成球体绕开了所有标定问题，法线真值已知，图像用同一个朗伯模型渲染，跑出来的误差基本只来自数值计算本身。源码包里的 scripts/make_sphere_dataset.py 就干这个事。
6.2 快速生成球体数据并评估误差
脚本生成一个半径 100 像素的球体，从 8 个不同方向打光，渲染出 8 张球体灰度图，同时保存真值法线和掩膜。生成后，用编译好的 photo_stereo 去处理这组数据，再把输出的 normal_map.exr 和真值法线比较，计算平均角度误差：
import cv2 import numpy as np gt = cv2.imread("sphere_gt_normal.exr", -1) est = cv2.imread("normal_map.exr", -1) mask = cv2.imread("sphere_mask.png", -1) > 0 dot = np.clip(np.sum(est * gt, axis=2), -1.0, 1.0) angle = np.degrees(np.arccos(dot)) valid = angle[mask] print("mean angle error:", np.mean(valid)) print("p95 angle error:", np.percentile(valid, 95))正常情况下平均角度误差应该在 1 度以内，p95 在 3 度以内。误差超过这个范围，优先检查光源方向文件的顺序和掩膜。合成球体因为是理想朗伯体，不存在高光和阴影干扰，一旦误差偏大，问题一定在代码层面。把这套验证跑通之后，再换真实物体数据，你就有了一个基线：新数据上的问题大概率是光源方向标定或输入图像质量造成的，而不是算法本身。从那以后，我每做一套光度立体数据，都会先强制走一遍合成球体验证，把代码层和物理层的误差分开看，这让我少踩了很多标定上的坑。希望帮到你。
本文还有配套的精品资源，点击获取
