光度立体三维重建——光源标定

一、光源标定问题

光源的标定是指 放置标定物等辅助物体估计光源的信息，包括光源的方向和强度等，光源信息的准确性极大影响重建的效果。光度立体的前提假设入射光为平行光，现实往往难以制造大面积平行光，所以通常的做法是点光源与物体的距离十倍于物体最大宽度时所发出的光可近似看成平行光。

Pentland最早提出计算光源方向，指出光源标定的线索更多地包含在图像的亮度梯度而不是图像的绝对亮度中。传统的标定中一般利用高光黑球来做标定物，通过点光源在黑球反射的高光信息与黑球的几何关系，求解光源方向。如图光源标定示意图，光源方向是 L，N 是球体的表面 法 向 量 ，视 角 方 向 为 V[0，0，1]，则 其 关 系 为 L = 2(N∙V)N - V ，V 已知，只要知道 N 即可求出 L，N 则可以通过高光黑球的已知半径 R 和光源照射在其表面的几何关系确定，r 可以通过图像处理的方法求得，z 可通过 r 和 R 求得。

光源标定的目的 是为了确定光源方向或强度，一般利用物体表面的明暗信息、阴影信息或者反射特性三类信息中的其中一种即可，也可以将三种信息结合。采用不同标定物，是利用不同标定物的反射特性获取更加精确的光源信息。光源标定方法分类总结如下图：

二、Matlab代码

close all;
% 定义目录和光源数量
directory = 'F:\......\data\4-901\';
numLights = 4; % 假设有 4 个光源

% 调用校准函数
L = CalibrateLights(directory, numLights);

function [L] = CalibrateLights(directory, numLights)
    
    % 检查目录路径是否以斜杠结尾，如果没有则添加斜杠。
    if directory(end) ~= '/' && directory(end) ~= '\'
        directory = [directory '/']; % 添加斜杠（支持 UNIX 和 Windows 风格路径）
    end

    % 初始化变量
    fileName = directory; % 初始化 fileName 为目录路径
    maskFileName = [fileName 'mask.bmp']; % 指定掩码文件的名称（修正了拼写错误）

    % 读取掩码图像
    try
        circle = imread(maskFileName); % 读取掩码图像
    catch ME
        error('无法读取掩码文件 %s，错误信息：%s', maskFileName, ME.message);
    end

    % 计算铬球的中心点
    maxval = max(circle(:)); % 找到图像中的最大值
    [circleRow, circleCol] = find(circle == maxval); % 找到最大值的所有坐标
    maxRow = max(circleRow); % 找到最大行坐标中的最大值
    minRow = min(circleRow); % 找到最小行坐标中的最小值
    maxCol = max(circleCol); % 找到最大列坐标中的最大值
    minCol = min(circleCol); % 找到最小列坐标中的最小值
    xc = double((maxCol + minCol)/2); % 计算列坐标的中心
    yc = double((maxRow + minRow)/2); % 计算行坐标的中心
    center = [xc, yc]; % 存储中心点坐标
    radius = double((maxRow - minRow)/2); % 计算半径（假设铬球是圆形的）

    % R: 反射方向，这里假设为垂直向上的方向。
    R = [1.0 1.0 1.0]; % 通常反射方向应为单位向量，例如 [0 0 1]

    % 初始化光源方向数组
    L = zeros(numLights, 3);

    % 计算光源方向
    for i = 1:numLights
        imgFileName = [fileName num2str(i-1) '.bmp']; % 根据索引构造图像文件名
        try
            image = imread(imgFileName); % 读取图像
        catch ME
            error('无法读取图像文件 %s，错误信息：%s', imgFileName, ME.message);
        end

        maxval = max(image(:)); % 找到图像中的最大值
        [pointRow, pointCol] = find(image == maxval); % 找到最大值的所有坐标
        nSize = size(pointRow, 1); % 计算最大值的数量
        px = sum(pointCol)/double(nSize); % 计算列坐标的平均值
        py = sum(pointRow)/double(nSize); % 计算行坐标的平均值
        Nx = px - xc; % 计算x方向上的差值
        Ny = -(py - yc); % 计算y方向上的差值（注意这里的负号）
        Nz = sqrt(radius^2 - Nx^2 - Ny^2); % 根据半径和Nx, Ny计算Nz
        if isnan(Nz) || Nz < 0
            error('计算出的Nz为负数或NaN，请检查输入数据。');
        end
        normal = [Nx, Ny, Nz]; % 计算法向量
        normal = normal / norm(normal); % 归一化法向量
        NR = normal(1)*R(1) + normal(2)*R(2) + normal(3)*R(3); % 计算法向量和反射方向的点积
        L(i,:) = 2*NR*normal - R; % 根据Lambertian反射模型计算光源方向
    end

    % 将新的光源方向写入测试文件
    fid = fopen('calibratedLight3250.txt', 'w'); % 打开文件以写入
    if fid == -1
        error('无法打开文件进行写入。');
    end
    fprintf(fid, '%d \n', numLights); % 写入光源数量
    for row = 1:numLights
        fprintf(fid, ' %10.5f %10.5f %10.5f \n', L(row,1), L(row,2), L(row,3)); % 写入光源方向
    end
    fclose(fid); % 关闭文件
end

三、实验效果

从四个方向逐次打光逐次采集图像，获得以下四张数据：

对以上数据进行数据融合并进行二值化操作（参考博客：OpenCV——图像融合）
获得以下数据：

输出标定矩阵：

参考：
[1]邓学良,何扬波,周建丰.基于光度立体的三维重建方法综述[J].现代计算机,2021,27(23):133-143.
[2]丁斌杰.面向工业场景的光度立体重建方法研究[D].浙江大学,2023.002432.
[3]基于光度立体视觉的三维重建方法
[4]Photometric Stereo 光度立体三维重建（四）——光源标定