# Interpolation by Solving a Minimization Problem (Optimization)

I will try to give the motivation behind this problem and later the math formality.
Given a grayscale image (1 Channel - $M \times N$ Matrix).
Someone marks some pixels as anchors.
Now, you need to interpolate the other pixels (Which are not anchors) by minimizing a given cost function s.t. the end result is an image which has the original image values at the anchors and interpolated values else were s.t. it minimizes the cost function.

Given an $M\times N$ matrix (A 1 channel image for that matter) $I$.

A subset of the elements (Pixels) in the matrix are marked as a reference and their location is a group marked as $S$.

The optimization cost function is given by:

$$E = \arg \min_{E} \sum_{\mathbf{r}} \left( E(\mathbf{r}) - \sum_{\mathbf{s} \in N(\mathbf{r})} {w}_{\mathbf{rs}} E(\mathbf{s}) \right)^2, \text{ s.t. } \forall p \in S \; E(\mathbf{p}) = I(\mathbf{p})\,,$$

where a bold letter $\mathbf{p}, \mathbf{r}, \mathbf{s}$ means an element (Pixel) location.

The group $N(r)$ is the neighborhood of $\mathbf{r}$, which is size $k$ namely, a $k$ by $k$ rectangle where $\mathbf{r}$ is in the middle.

The weights ${w}_{\mathbf{rs}}$ are defined as following:

$${w}_{\mathbf{rs}} \propto \exp\left(-\frac{(I(\mathbf{r}) - I(\mathbf{s}))^2}{2\sigma^2_r} \right )\ \ \text{ s.t. } \sum_{\mathbf{s} \in N(\mathbf{r})} w_{\mathbf{rs}} = 1\,.$$

Namely, the weights are normalized to 1 within the neighborhood. The variance is calculated on the matrix $I$ in the neighborhood (You can assume it is given).

So the problem is to find a matrix $E$ which is equal to $I$ on all reference points and interpolates other places by bringing the cost function to the minimum.

It looks like a weighted least squares per neighborhood, yet I couldn't formalize it in a way that can be easily calculated and solved in e.g. MATLAB.

• What is your difficulty? What have you tried already and what do you mean by "I couldn't formalize it a way that can be easily calculated and solved"? If you write down the optimality conditions, they can be written as solving a linear system with a matrix that has $k^2$ entries per row. Where in this process are you encountering problems? Apr 15 '14 at 0:08
• @WolfgangBangerth, Could you assist me even with the ${k}^{2}$ entries solution. Just the how to formalize it? Later maybe we could think about faster and better. I'm not experienced with those kind of things. Thank You.
– Royi
Apr 15 '14 at 0:11
• @WolfgangBangerth, I meant in the typical WLS problem I know we have the measurements and the residuals are well defined. Here I have no idea how define those. Let's say I built a matrix in each row the corresponding weights, as you wrote, ${k}^{2}$, Now it should be multiplied by the model parameters, assuming those are the pixel, then what's the result? It can be the image, it also can't be the anchors (Size doesn't match). I guess I miss it.
– Royi
Apr 15 '14 at 1:41
• @WolfgangBangerth, Something doesn't add up with ${k}^{2}$ elements in a row. Since if you solve the linear equation of such a matrix you get the solution in the neighborhood only. While it has to be a global (N by M).
– Royi
Apr 15 '14 at 8:17
• You have it all there already. You have an optimization problem with constraints. Write down the Lagrangian as always, take the derivatives to form the optimality conditions, and you have a linear system for all of your unknowns. Apr 15 '14 at 12:30

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;

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;

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:

1. 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.
2. 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$.
3. 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.
4. 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.

• Tell me, you said it is a little different optimization problem, could you say what exactly? Can the original be solved as well?
– Royi
Apr 16 '14 at 7:14
• By the way, at the page you linked they have a "Solver". Though I don't understand what's going on there.
– Royi
Apr 16 '14 at 8:19
• @Drazick updated my answer to include practical MATLAB code. This code can be improved but hopefully it captures the main ideas. Apr 19 '14 at 7:11
• Tried your code and it works! Got 2 things: 1. Look down at the derivation I made into Quadratic form of the original, Do you think you can create an optimized code for it? Could you say something about solving it for large images under memory constraints? 2. Assuming we have 2 optimized codes, one using you problem the other for mine, which should be faster? Namely, Is any of the forms inherently easier to solve?
– Royi
Apr 21 '14 at 22:30

Using few tricks the equation can be transformed into classic Quadratic Programming problem:

$$E = \arg \max_{E} \sum_{\mathbf{r}} \left( E(\mathbf{r}) - \sum_{\mathbf{s} \in N(\mathbf{r})} {w}_{\mathbf{rs}} E(\mathbf{s}) \right)^2$$

Defining $\boldsymbol{e} = vec \left ( E \right )$, Namely vectorizing the matrix, and defining a matrix $A$ defined the proper usage of the weights ${w}_{\mathbf{rs}}$, Namely put them in the correct entry at each row per pixel the above could be rewritten as:

$$\mathbf{e} = \arg \max_{\mathbf{e}} {\left ( \mathbf{e} - A\mathbf{e} \right ) }^{T} \left ( \mathbf{e} - A\mathbf{e} \right ) = \arg \max_{\mathbf{e}} {\left( \left ( I - A \right ) \mathbf{e} \right )}^{T} \left( \left ( I - A \right ) \mathbf{e} \right )$$

Defining $B = I - A$ yields:

$$\mathbf{e} = \arg \max_{\mathbf{e}} {\left( \left ( I - A \right ) \mathbf{e} \right )}^{T} \left( \left ( I - A \right ) \mathbf{e} \right ) = \arg \max_{\mathbf{e}} {\mathbf{e}}^{T} {B}^{T} B \mathbf{e}$$

Defining $L = {B}^{T} B$ , namely $L$, is a Symmetric Positive Definite Matrix yields:

$$\mathbf{e} = \arg \max_{\mathbf{e}} {\mathbf{e}}^{T} L \mathbf{e}$$

Adding the equality constraints using $C \mathbf{e} = \mathbf{t}$ where $\mathbf{t}$ is vectorization of all "Reference Pixels"

All in all this yields a Classic Quadratic Programming problem given by:

$$\mathbf{e} = \arg \max_{\mathbf{e}} {\mathbf{e}}^{T} L \mathbf{e}, \; s.t. \; C \mathbf{e} = \mathbf{t}$$

Few notes:

1. The matrix $L$ is Symmetric and Positive Definite matrix and usually sparse (Assuming "Small Neighborhood").
2. The matrix $L$ wast built by the Matrix $I - A$ which each row of it was normalized to be zero (Each row of $A$ is normalized to 1 ).

Now two things are needed:

1. Optimized fast algorithm to solve this problem utilizing its properties.
2. Assuming the image is large, are there methods to solve this problem by blocks to avoid memory limits?