
简介这是一份基于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 或直接保存的 TIFFJPG 压缩在高光和暗部容易引入噪点影响法线精度。光源方向文件的格式很简单每行三个浮点数分别代表 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_TYPERelease .. 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::vectorcv::Mat readImages(const std::vectorstd::string paths) { std::vectorcv::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 一次求出所有像素的 GG 是 3 行 M 列每一列就是该像素的 ρ·n。void solveNormalMap(const std::vectorcv::Mat images, const std::vectorcv::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.atdouble(i, 0) lights[i][0]; L.atdouble(i, 1) lights[i][1]; L.atdouble(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.atuchar(r, c) 0) { normalMap.atcv::Vec3f(r, c) cv::Vec3f(0.0f, 0.0f, 1.0f); albedoMap.atdouble(r, c) 0.0; continue; } double gx G.atdouble(0, idx); double gy G.atdouble(1, idx); double gz G.atdouble(2, idx); double len std::sqrt(gx * gx gy * gy gz * gz); albedoMap.atdouble(r, c) len; if (len 1e-12) { normalMap.atcv::Vec3f(r, c) cv::Vec3f( (float)(gx / len), (float)(gy / len), (float)(gz / len)); } else { normalMap.atcv::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.atcv::Vec3f(r, c); cv::Vec3b out; out[0] cv::saturate_castuchar((n[2] 1.0) * 127.5); out[1] cv::saturate_castuchar((n[1] 1.0) * 127.5); out[2] cv::saturate_castuchar((n[0] 1.0) * 127.5); vis.atcv::Vec3b(r, c) out; } } return vis; }OpenCV 里 Mat 的通道顺序是 BGR所以这里 out[0] 对应蓝色通道保存的是法线 z 分量。这样做的结果是正对相机的平面法线是 (0,0,1)映射后蓝色通道是 255看起来就是一片蓝色这是法线贴图的标准观感。out[2] 对应红色通道保存的是 x 分量。如果你要用 Blender、Unity 这类工具读取法线贴图它们内部默认的通道约定可能是 Rx、Gy、Bz那就要做一次通道重排。源码里把原始浮点法线也保存了一份避免 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 字节约 64MB8 张就是 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, axis2), -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 度以内。误差超过这个范围优先检查光源方向文件的顺序和掩膜。合成球体因为是理想朗伯体不存在高光和阴影干扰一旦误差偏大问题一定在代码层面。把这套验证跑通之后再换真实物体数据你就有了一个基线新数据上的问题大概率是光源方向标定或输入图像质量造成的而不是算法本身。从那以后我每做一套光度立体数据都会先强制走一遍合成球体验证把代码层和物理层的误差分开看这让我少踩了很多标定上的坑。希望帮到你。本文还有配套的精品资源点击获取