This is the classic colorization using optimization problem.
Optimization and Linear Algebra
To see how this can be expressed as a linear system, it's helpful to use a slightly different notation (and a slightly different objective function). Think of your image as graph $G$ with a node for each pixel in the image. There is an edge $(i,j)$ between two pixels $i$ and $j$, if these pixels are neighbors in the image.
Let $x_i$ be the intensity of pixel $i$, then you can write the problem as
$
\begin{array}{ll}
\underset{x}{\mathrm{minimize}} & \frac{1}{2} \sum_{(i,j) \in E} w_{ij} (x_i - x_j)^2 \\
\mathrm{subject\ to} & x_i = d_i, i \in S.
\end{array}
$
For each edge in the graph we have weights denoted by $w_{ij} = w_{ji}$. Given this weighted graph $G(V, E, w)$ we can define the adjacency matrix $A$ by
$ A_{ij} = \begin{cases} w_{ij} & \text{if}\ (i,j) \in E\\ 0 & \text{otherwise} \end{cases} $
The degree matrix of a weighted graph is a diagonal matrix $D$ such that
$ D_{ii} = \sum_{j} A_{ij}$
The Laplacian matrix of a weighted graph is defined as $L = D - A$. With this definition of the Laplacian we can see that
$ x^T L x = \sum_{(i,j) \in E} w_{ij} (x_i - x_j)^2$
Thus we can rewrite the above problem in matrix form as:
$\begin{array}{ll}
\underset{x}{\mathrm{minimize}} & \frac{1}{2} x^T L x \\
\mathrm{subject\ to} & x_m = d\\
\end{array}
$
Here $L$ is the Laplacian matrix for the image, $x_m$ is the set of unknowns corresponding to marked pixels, and $d$ is the vector of annotated values (anchors).
Let $P$ be a permutation matrix that orders $x$ so that the marked and unmarked unknowns are grouped such that
$ Px = \begin{bmatrix} x_m \\ x_u \end{bmatrix}$
Then we have that
$ \frac{1}{2} x^T L x = \frac{1}{2} x^T P^T (PLP^T) Px = \frac{1}{2} \begin{bmatrix} x_m^T & x_u^T \end{bmatrix} \begin{bmatrix} L_m & R^T \\ R & L_u \end{bmatrix} \begin{bmatrix} x_m \\ x_u \end{bmatrix}$
Here we have that $L_u = L_u^T \succ 0$. Since the values of the variables $x_m$ are fixed ($x_m = d$) we need only optimize over the variables $x_u$. Thus we can solve the unconstrained problem
$ \underset{x_u}{\mathrm{minimize}}\ F(x_u) = \frac{1}{2} x_u^T L_u x_u + x_m^T R^T x_u$
Since $L_u \succ 0$, the unique minimizer $x_u^\star$ is characterized by $\nabla F(x_u^\star) = 0$. We have that $\nabla F(x_u) = L_u x_u + R x_m$. Thus we need only solve the linear symmetric positive-definite system
$L_u x_u = -Rd$
You can solve this with a direct method (i.e. \ in MATLAB or conjugate gradient). For a class in grad school I implemented PCG in CUDA on a graphics card to color a 1050 x 3360 image which corresponded to a linear systems with 2.1 million unknowns. It took 230s on two Quadro FX 5600 GPUs to solve the system.
Sample MATLAB code
Below is sample MATLAB code which computes a colorization using the above images
original = double(imread('example.bmp'))/255;
marked = double(imread('example_marked.bmp'))/255;
annotate_thresh = 0.25; % may need to be adjusted
colorim = sum(abs(original - marked), 3) > annotate_thresh;
colorim = double(colorim);
spy(colorim);
anchorind = find(colorim);
[m,n,p] = size(marked);
assert(p == 3);
N = m*n;
% convert from RGB to NTSC using the matrix T
T = [1.0 0.956 0.621; 1.0 -0.272 -0.647; 1.0 -1.106 1.703];
yiq = reshape(marked(:), N, 3)/T';
yiq = reshape(yiq, m, n, 3);
chromi = yiq(:,:,2);
chromq = yiq(:,:,3);
% Get grayscale value from original image
yiq = reshape(original(:), N, 3)/T';
yiq = reshape(yiq, m, n, 3);
gray = yiq(:,:,1);
% Get anchors from the I and Q color components
anchor(:,1) = chromi(anchorind);
anchor(:,2) = chromq(anchorind);
edges=[(1:N)' ((1:N)+1)']; %add downard edge between node i with node i+1
edges=[edges; (1:N)' (1:N)'+m]; %add leftward edge between node i and i+m
% exclude edges which link to nodes outside the domain
excluded=edges(:,1)>N | edges(:,1)<1 | edges(:,2)>N |edges(:,2)<1;
edges(excluded, :) = [];
% remove edge from bottom of column to top of next column
edges((m:m:(m-1)*n)',:)=[];
beta = 200;
EPSILON = 1e-5;
% calculate edge weights
graydistance = abs(gray(edges(:,1)) - gray(edges(:,2)));
ming = min(graydistance); % it's helpful to normalize distances
maxg = max(graydistance);
graydistance = (graydistance - ming)./(maxg - ming);
weights=exp(-(beta*graydistance)) + EPSILON;
%Build sparse weighted adjacency matrix
W=sparse([edges(:,1);edges(:,2)],[edges(:,2);edges(:,1)], ...
[weights;weights],N,N);
% Build Degree matrix
D = diag(sum(W));
% Build Laplacian
L=D-W;
mark = anchorind;
rest = 1:n*m;
rest(mark) = [];
rest = rest';
% Extract the portion of the Lapalcian matrices
Lu = L(rest, rest);
R = L(rest, mark);
recover = zeros(m, n, 3);
recover(:,:,1) = gray;
xchannel = zeros(m*n,1);
for i=1:2
d = anchor(:,i);
xrest = -Lu\(R*d);
xchannel(rest) = xrest;
xchannel(mark) = d;
recover(:,:,i+1) = reshape(xchannel, m, n);
end
% Convert from YIQ back to RGB using T matrix
recover = reshape(recover, m*n, 3);
recover = recover*T';
% Keep RGB values within [0,1]
recover = min(max(0, recover), 1);
final = reshape(recover, m, n, 3);
figure(1);
image(final);
The output is the image below
Why a different objective function
I choose a different objective function than the one stated in your problem.
There are several reasons for this:
- The quantity $w_{ij} (x_i - x_j)^2$ has a clear interpretation in the objective. You pay a weighted penalty based on the square of the difference
between two adjacent pixels.
- Using the definition of the Laplacian matrix $L$ it's easy to go from this objective function to the matrix equation $x^T L x$.
- Using the definition of $L = A^T C A$ (where $A$ is the edge-node adjacency matrix and $C$ is a diagonal matrix containing the edge weights) it's easy to see that the Laplacian is positive semi-definite. If we add a small bit of regularization to make it positive definite, we can solve systems with $L$ via a Cholesky factorization or preconditioned conjugate gradient.
- It's possible to do everything with the objective function you defined above. It's just more difficult to figure out how to translate it into a matrix equation.
Note the idea of using this different objective function is not my own. I found it in a paper by Leo Grady on Random Walks for Image Segmentation (See presentation Random Walks for Image Segmentation).
Working with large images
An image with $m$ rows and $n$ columns will yield a matrix $L$ of size $mn \times mn$. The number of nonzeros in $L$ is related to the size of the neighborhood around a pixel. In the code above, each pixel can have up to four neighbors: up, down, left, right. So the number of nonzeros in $L$ is bounded by $5mn$. Holding $L$ in memory is equivalent to having a small number of copies of the image in memory.
If you want to work with large images, you probably want to use some sort of iterative method. Iterative methods don't require you to form $L$ directly, but rather compute matrix-vector products. That is given a vector $x$ compute a vector $y$ such that $y=Lx$. As I mentioned above, I used preconditioned conjugate gradient, with diagonal preconditioning, on multiple graphics cards to solve this problem. This worked OK, but was still fairly slow on large images. To work with large images efficiently you may need to come up with a better preconditioner. There's been some research on this, for example see Richard Szeliski's paper on "Locally Adapted Hierarchical Basis Preconditioning". I think others have also investigated using multigrid methods to solve this problem.
