计算机视觉：朗伯光度立体法（Lambertian Photometric Stereo）
光度立体法简介
光度立体法，即Photometric Stereo, 最早是由当时在MIT的人工智能实验室的Robert J. Woodham教授在1978年左右提出。他在1979年的论文《Photometric stereo: A reflectance map technique for determining surface orientation from image intensity》，以及1980年的论文《Photometric Method for Determining Surface Orientation from Multiple Images》中比较系统的阐述了整套理论框架。光度立体法可以根据二维纹理信息提取出三维模型，其典型应用是检测物体表面微小变化。
朗伯光度立体法基于Woodham算法，Woodham在论文中提出三个基本假设。在这三个基本假设下完成整套理论框架的推演。这三个基本假设分别是：
- 假定相机是无畸变成像，也就是说必须使用远心镜头或者长焦镜头。
- 假定每一个光源发射的光束都是平行且均匀的，也就是说必须使用具有均匀强度的远心照明光源，或者使用远距离的点光源代替。
- 物体必须具有朗伯（lambertian）反射特性，即它必须以漫反射的方式反射入射光。
朗伯光度立体法的大致思路是：
当相机和目标物体相对位置固定不变时，使用不同方向的光源照射同一目标物体，相机可以拍摄到目标物体带有不同明暗分布的图像（至少需要三张图），再通过求解基于朗伯反射原理的反射方程组，求解目标表面的法向分布或者albedo图。
带有远心镜头的相机必须与被测物体表面垂直安装，在采集多幅图像时，一定要保证相机和物体不被移动。相反，对于采集至少三张的灰度图像，其每次取像的照明方向必须改变（相对于相机）。对于采集的多张图像中的每一幅图，照明方向必须指定Slants和Tilts两个参数角度，其描述了相对于当前场景的光照角度。为了更好的理解这两个参数含义，我们假定光源射出的光束是平行光，镜头是远心镜头，相机垂直于物体表面。
朗伯光度立体法算法原理
根据Lambert模型：
I
=
ρ
L
⋅
N
\textbf{I} = \rho\ \textbf{L}\cdot\textbf{N}
I=ρ L⋅N
式中，
ρ
\rho
ρ 为表面反射率（albedo），其值介于0 - 1之间，反映了物体表面特性；
N
\textbf{N}
N 为表面法线（normal）；
L
\textbf{L}
L 为照明方向。将
ρ
\rho
ρ 和
N
\textbf{N}
N 合并起来用
G
\textbf{G}
G来表示，则有：
I
=
L
T
G
.
\textbf{I} = \textbf{L}^T\textbf{G}.
I=LTG.
每个像素位置对应的
G
\textbf{G}
G是三维向量（方向为normal，模长为albedo），因此单个光源无法求解该方程，至少需要三个不同位置的光源，求解得到三个未知量。
由最小二乘法：
min
G
∣
∣
I
−
L
T
G
∣
∣
2
\mathop{\min}\limits_{\textbf{G}} ||\textbf{I}-\textbf{L}^T\textbf{G}||^2
Gmin∣∣I−LTG∣∣2
可以求得：
G
=
(
L
L
T
)
−
1
LI
\textbf{G} = (\textbf{L}\textbf{L}^T)^{-1}\textbf{L}\textbf{I}
G=(LLT)−1LI
我们求得
G
\textbf{G}
G 的模长就是albedo:
ρ
=
∣
∣
G
∣
∣
\rho = ||\textbf{G}||
ρ=∣∣G∣∣
G
\textbf{G}
G 归一化后的单位向量矩阵就是normal:
N
=
G
∣
∣
G
∣
∣
\textbf{N} = \frac{\textbf{G}}{||\textbf{G}||}
N=∣∣G∣∣G
朗伯光度立体法matlab程序
注：此代码只涉及核心算法，不包含数据的读入，与结果的输出。
function[Albedo, Normal, Re_rendered]=Photometric_Stereo(data)assert(size(data.imgs,1)==size(data.s,1),'Size mismatched!');
num =size(data.imgs,1);% Get image dimensions
im = data.imgs{1};[im_h, im_w,~]=size(im);% Initialize T, a im_h-by-im_w-by-num matrix, whose (h, w,:) holds
% the intensities at (h, w)for all p different lightings
im_T =zeros(im_h, im_w, num);% Initialize im_R, a im_h-by-im_w-by-num matrix
im_R =zeros(im_h, im_w, num);% Initialize im_G, a im_h-by-im_w-by-num matrix
im_G =zeros(im_h, im_w, num);% Initialize im_B, a im_h-by-im_w-by-num matrix
im_B =zeros(im_h, im_w, num);% For each image
for idx =1:num
im = data.imgs{idx};
imGray =rgb2gray(im);% Loop thru each pixel
for h =1:im_h
for w =1:im_w
% If in the mask
if data.mask(h, w)im_R(h, w, idx)=im(h, w,1);im_G(h, w, idx)=im(h, w,2);im_B(h, w, idx)=im(h, w,3);im_T(h, w, idx)=imGray(h,w);
end
end
end
end
% Initialize Normal, a im_h-by-im_w-by-3 matrix
Normal =zeros(im_h, im_w,3);% Initialize Albedo, a im_h-by-im_w-by-1 matrix
Albedo =zeros(im_h, im_w,3);% Initialize Re_rendered, a im_h-by-im_w-by-3 matrix
Re_rendered =zeros(im_h, im_w,3);% Loop thru each location
for h =1:im_h
for w =1:im_w
% If in the mask
if data.mask(h, w)%% Normal
% Lightings
L = data.s;% Intensities
I =reshape(im_T(h, w,:),[num,1]);% Dealing with shadows and highlights
for k =1:10
max_index =find(I==max(I));I(max_index)=[];L(max_index,:)=[];
min_index =find(I==min(I));I(min_index)=[];L(min_index,:)=[];
end
% Solve surface normals and albedo
G =(L.'*L)\(L.'*I);ifnorm(G)~=0% Normalize n
n = G./norm(G);else
n =[0;0;0];
end
% Save
Normal(h, w,:)= n;%% Albedo
% a_R
L = data.s;
I_R =reshape(im_R(h, w,:),[num,1]);for k =1:10
max_index =find(I_R==max(I_R));I_R(max_index)=[];L(max_index,:)=[];
min_index =find(I_R==min(I_R));I_R(min_index)=[];L(min_index,:)=[];
end
G_R =(L.'*L)\(L.'*I_R);
a_R =norm(G_R);Albedo(h, w,1)= a_R;% a_G
L = data.s;
I_G =reshape(im_G(h, w,:),[num,1]);for k =1:10
max_index =find(I_G==max(I_G));I_G(max_index)=[];L(max_index,:)=[];
min_index =find(I_G==min(I_G));I_G(min_index)=[];L(min_index,:)=[];
end
G_G =(L.'*L)\(L.'*I_G);
a_G =norm(G_G);Albedo(h, w,2)= a_G;% a_B
L = data.s;
I_B =reshape(im_B(h, w,:),[num,1]);for k =1:10
max_index =find(I_B==max(I_B));I_B(max_index)=[];L(max_index,:)=[];
min_index =find(I_B==min(I_B));I_B(min_index)=[];L(min_index,:)=[];
end
G_B =(L.'*L)\(L.'*I_B);
a_B =norm(G_B);Albedo(h, w,3)= a_B;%% Re_rendered
Re_rendered(h, w,:)=[a_R;a_G;a_B]*dot(n,[0;0;1]);
end
end
end
完整程序请移步至此链接下载
示例
Albedo图
反照率图：
Normal图
以RGB线性编码的法线贴图：
Re_rendered图
在照明方向和观察方向一致时，利用normal和albedo图重新渲染的图片：
参考文献
光度立体简介
Halcon 光度立体法（photometric_stereo）详解z
版权归原作者 cv-player 所有， 如有侵权，请联系我们删除。