How do you formulate the linear least-squares method for radiometric calibration?

Question

In Debevec and Malik (mentioned similarly in Forsyth and Ponce's Computer Vision: A Modern Approach) they highlight a method of solving the camera response function using linear least-squares.

We collect image intensity data for a number of points with

$$ I_{pk} = f(E_{p}\Delta t_{k}) $$

Where $I_{pk}$ is the image intensity for pixel $p$ for the $k$-th exposure time
$t_{k}$ is the $k$-th exposure time
$E_{p}$ is the intensity of the surface projected onto the image pixel
$f$ is the camera response function

$I_{pk}$ and $t_{k}$s are known and we are using least-squares to solve for the unknown $E_{p}$'s and $f$. This problem is turned linear by taking $g = \ln f^{-1}$ and taking the logs we can formulate the above as:

$$g(I_{pk}) = \ln(E_{p}) + ln(\Delta t_{k}) $$

Then minimizing this as the following objective function by choice of $g$:

$$\sum(g(I_{pk}) - \ln(E_{p}) + ln(\Delta t_{k}))^2 + \sum_{z} g''(z)^2$$

Where $z$ is the discrete domain of all pixel intensities, and $g''(z) = g(z - 1) - 2g(z) +g(z + 1)$ creates a smoothness penalty on $g$.

From Debevec and Malik

Because it is quadratic in the $E_{i}$’s and $g(I)$’s, minimizing [the above] is a straightforward linear least squares problem

I don't understand how we've reached the following conclusion. As mentioned by the paper, we should be able to solve the following by creating the correct matrix and performing singular-value decomposition. What is the construction of this matrix?

Are we saying that $g = \ln f^{-1}$ is a linear function?

Answer 1

This is a Penalized Least Squares problem, with penalty on the second derivative of $g()$ (so we are not saying that the function is linear). In practice you have one unknown parameter for each pixel intensity $I_{pk}$, that represents the value of the function at that point $g(I_{pk})$. You just have to express your penalty in a matrix form $g(I)^T M g(I)$, where $M$ is a tridiagonal matrix with -2s on the main diagonal, and 1s above and below the main diagonal (looks right, but not completely sure about it). Then either just give it to an optimizer to estimate the values $g(I)$ and $E$ at the same time or you estimate them one at the time. For example for fixed $E$, the closed form solution for $g(I)$ should be $\hat{g}(I)=(I^TI + M)^{-1}I^T(ln(\Delta t) -ln(E))$. Then for $E$ you can just use an optimizer (with fixed $\hat{g}(I)$) and alternate until convergence.