How to understand the implementation codes in Debevec and Malik's HDR paper

Question

The paper titled "Reconvering High Dynamic Range Radiance Maps from Photographs" from Debevec and Malik is a seminal paper in high dynamic imaging. In this paper, it proposed a method to estimate the camera response function (CRF), which is equal to optimize the following objective function:

$$\mathcal O = \sum_{i=1}^N\sum_{j=1}^P\left\{w\left(Z_{ij}\right)\left[g\left(Z_{ij}\right)-\ln E_i-\ln\Delta t_j\right]\right\}^2+\lambda \sum_{z=Z_{\rm min}+1}^{Z_{\rm max}-1}\left[w(z)g''(z)\right]^2$$

The author also provides MATLAB codes to optimize this object function:

%
% gsolve.m − Solve for imaging system response function
%
% Given a set of pixel values observed for several pixels in several
% images with different exposure times, this function returns the
% imaging system’s 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
% w(z) is the weighting function value for pixel value z
%
% Returns:
%
% g(z) is the log exposure corresponding to pixel value z
% lE(i) is the log film irradiance at pixel location i
%
function [g,lE]=gsolve(Z,B,l,w)

n = 256;

A = zeros(size(Z,1)*size(Z,2)+n+1,n+size(Z,1));
b = zeros(size(A,1),1);

%% Include the data−fitting equations

k = 1;
for i=1:size(Z,1)
  for j=1:size(Z,2)
    wij = w(Z(i,j)+1);
    A(k,Z(i,j)+1) = wij; A(k,n+i) = −wij; b(k,1) = wij * B(i,j);
    k=k+1;
  end
end

%% Fix the curve by setting its middle value to 0

A(k,129) = 1;
k=k+1;

%% Include the smoothness equations

for i=1:n−2
  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

%% Solve the system using SVD

x = A\b;

g = x(1:n);
lE = x(n+1:size(x,1));

The implementation codes are well understood, and my question is why this implementation can optimize the object function. If I understand well about this implementation, the solution to the object function is the one that makes the gradients of the object function equal to zero. However, one thing I am confused is that why the first part and the second part of the object function are treated independently. How could it be possible?

Answer 1

OK, this is incomplete but this is how far I've gotten:

My previous answer, which you can see in the edits was incorrect.

Your question boils down to:

However, one thing I am confused is that why the first part and the second part of the object function are treated independently.

So let's try having a look at the objective function: $$ \begin{array}{rcl} {\cal O} &=& \displaystyle \sum_{i=1}^{N} \sum_{j=1}^P \left\{ w(Z_{ij}) [ g(Z_{ij}) - \ln E_i - \ln \Delta t_j ] \right\}^2 + \lambda \sum_{z =Z_{\rm min}+1}^{Z_{\rm max}} [w(z) g^"(z)]^2\\ &=& \displaystyle \sum_{i=1}^{N} \sum_{j=1}^P w(Z_{ij})^2 [ g(Z_{ij})^2 + (\ln E_i)^2 + (\ln \Delta t_j)^2 - 2 g(Z_{ij}) \ln E_i - 2 g(Z_{ij}) \ln \Delta t_j + 2 \ln E_i \ln \Delta t_j ]\\ &&+ \lambda \displaystyle \sum_{z =Z_{\rm min}+1}^{Z_{\rm max}} [w(z) \{ g(z-1) - 2g(z) + g(z + 1) \} ]^2 \end{array} $$ where the second line just expands the first term and substitutes for $g^"$ as the paper suggests.

Now, we need to choose the $\ln E_i$ and $g$ values to minimize this so we must take the derivative of $\cal O$ with respect to $\ln E_i$ and to $g$.

$$ \frac{d \cal O}{d \ln E_i} = \sum_{i=1}^{N} \sum_{j=1}^P w(Z_{ij})^2 [ 2 \ln E_i - 2 g(Z_{ij}) + 2 \Delta t_j ] $$