How to understand the implementation codes in Debevec and Malik's HDR paper 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?
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} &amp;=&amp; \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\\
&amp;=&amp; \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
 ]\\
&amp;&amp;+ \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
 ]
$$