相机响应曲线计算（Debeve and Malik）

1.背景

相机的灰度值与实际辐照度值的关系不成线性关系，因此需要进行预标定。

2.算法（Debeve and Malik）

图片引自知乎文章，横坐标表示曝光度ln（E*t），因为辐照度E不知道，所以假设为1，纵坐标表示像素值灰度百分比，曲线表示一个像素在不同曝光时间下的变化。上述曲线假设辐照度E都为1，所有这几条曲线没有重合，但理论上这几条曲线应该重合在一起，因此如果获取到这些像素的真实辐照度，即将这些曲线左右平移,表示ln(Et)=ln(E)+LN(1t)，重合在一起，即可得到真实照度与像素灰度值关系。因为曲线重合在一起后会有多个结果，为了消除全局偏移不确定性，可以将中间强度的像素值的曝光度置为0进行确定。

上述移动的本质其实就是求解实际辐照度，具体推导思路如下所示。

其中Z表示像素灰度值，i表示像素序号，j表示曝光时间序号，$Z_{ij}$表示曝光时间为$△t$时第i个像素的灰度值；E表示辐照度；f表示转换函数。可以推导为



g表示转换函数，将像素强度转变为曝光度ln（E*t）。

通过上述公式，可以知道要想总曲线达到最佳重合，可以列出公式

其中O为需要最小化的能量函数。前面一项为函数的最小二乘误差，后一项表示曲线的光滑度，越小表示曲线越光滑，其中$\lambda$表示光滑因子。由于像素灰度值在中间范围时较为可靠，因此可以考虑添加权重函数$\omega$，表示为

则能量函数可转换为

这一公式中需要求解的未知数为256+N个（256表示g函数的一一对应，N表示N个像素点的辐照度），可以列出的方程为NP个，同时由于上述将中间灰度值的辐照度置为0，则只要满足256+P$\le$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));