1.背景
相机的灰度值与实际辐照度值的关系不成线性关系,因此需要进行预标定。
2.算法(Debeve and Malik)
图片引自知乎文章,横坐标表示曝光度ln(E*t),因为辐照度E不知道,所以假设为1,纵坐标表示像素值灰度百分比,曲线表示一个像素在不同曝光时间下的变化。上述曲线假设辐照度E都为1,所有这几条曲线没有重合,但理论上这几条曲线应该重合在一起,因此如果获取到这些像素的真实辐照度,即将这些曲线左右平移,表示ln(Et)=ln(E)+LN(1t),重合在一起,即可得到真实照度与像素灰度值关系。因为曲线重合在一起后会有多个结果,为了消除全局偏移不确定性,可以将中间强度的像素值的曝光度置为0进行确定。
上述移动的本质其实就是求解实际辐照度,具体推导思路如下所示。
其中Z表示像素灰度值,i表示像素序号,j表示曝光时间序号,
Z
i
j
Z_{ij}
Zij表示曝光时间为
△
t
\triangle t
△t时第i个像素的灰度值;E表示辐照度;f表示转换函数。可以推导为
g表示转换函数,将像素强度转变为曝光度ln(E*t)。
通过上述公式,可以知道要想总曲线达到最佳重合,可以列出公式
其中O为需要最小化的能量函数。前面一项为函数的最小二乘误差,后一项表示曲线的光滑度,越小表示曲线越光滑,其中
λ
\lambda
λ表示光滑因子。由于像素灰度值在中间范围时较为可靠,因此可以考虑添加权重函数
ω
\omega
ω,表示为
则能量函数可转换为
这一公式中需要求解的未知数为256+N个(256表示g函数的一一对应,N表示N个像素点的辐照度),可以列出的方程为NP个,同时由于上述将中间灰度值的辐照度置为0,则只要满足256+P
≤
\leq
≤NP+1即可进行解算。
3.matlab代码
相关求解代码如下
function [g,lE]=gsolve(Z,B,l)
% Given a set of pixel values observed for several pixels in several
% images with different exposure times, this function returns the
% imaging system誷 response function g as well as the log film irradiance
% values for the observed pixels.
%
% Assumes:
%
% Zmin = 0
% Zmax = 255
%
% Arguments:
%
% Z(i,j) is the pixel values of pixel location number i in image j
% B(j) is the log delta t, or log shutter speed, for image j
% l is lamdba, the constant that determines the amount of smoothness
%
% Returns:
%
% g(z) is the log exposure corresponding to pixel value z ,g=log(E)+log(t)
% lE(i) is the log film irradiance at pixel location i
% w(z) is the weighting function value for pixel value z
% here we set w(z) to a constant
n = 256;
m = size(Z,1);
p = size(Z,2);
w = ones(n,1)/n;
% 创建稀疏零矩阵
A = sparse([], [], [], m*p+n+1,n+m, m*p*2 + n*3);% 前面两个变量为
b = zeros(size(A,1),1);
k = 1; %% Include the data-fitting equations
for i=1:m
for j=1:p
wij = w(Z(i,j)+1);
A(k,Z(i,j)+1) = wij;
A(k,n+i) = -wij;
b(k,1) = wij * B(j);
k=k+1;
end
end
A(k,129) = 1; %% Fix the curve by setting its middle value to 0
k=k+1;
% 计算g(Z_ij)前后的梯度
for i=1:n-2 %% Include the smoothness equations
A(k,i)=l*w(i+1);
A(k,i+1)=-2*l*w(i+1);
A(k,i+2)=l*w(i+1);
k=k+1;
end
% 相当于Ax=b,其中x前256为为g(i,j),后面为E_i,即对应像素的辐照度。
% A前面 m*p行每行有两个元素,分别是wij和-wij,
% 后面1行是129列处值为1,最后的255行为三个数w,-2w,w,表示梯度计算
%
% b 为wij * B(j)
x = A\b; %% Solve the system using SVD
g = x(1:n);
lE = x(n+1:size(x,1));