# Interpolation Optimization Problem

I will try to give the motivation behind this problem and later the math formality.
Given a grayscale image (1 Channel - M by 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$ by $N$ matrix (A 1 channel image for that matter) $I$.

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

The optimization cost function is given by:

$$\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 minimum.

It looks like a weighted least squares per neighborhood (The inner brackets).

I couldn't formalize it (For the whole matrix) a way that can be easily calculated and solved in e.g. MATLAB. Is there a way to formalize it as classic Weighted LS problem?
Or any other form which using the classic tools will bring solution (Get the interpolated matrix)?

Could any one help with that?

P.S. I tried using the Weighted LS per neighborhood disregarding the rest didn't yield the expected results (Obviously).
Tried it just to see how far the real solution is from this naive solution.

• Have you tried using a QP solver like CPLEX? It seems you could restate this as a min-cost flow problem. – Nicholas Mancuso Apr 15 '14 at 2:16
• Could you elaborate on that? I am not familiar with it. – Royi Apr 15 '14 at 2:18
• You could state this as a min-cost flow, where the cost over an edge is the pseudo-guassian you have. Instead of conserved flow you are trying to minimize the squared difference between in-flow and out-flow. – Nicholas Mancuso Apr 15 '14 at 2:43
• What do you mean by "expected results"; can you show an example? Are you sure the implementation and not the formalization is the issue? – Raphael Apr 15 '14 at 10:22
• I meant that the solution isn't by solving the inner brackets as classic WLS problem. Since the same element "plays" different role when $\mathbf{r}$ changes. The problem is how to formalize this function in a matrix form which can be minimized. – Royi Apr 15 '14 at 11:08

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 $\mathbf{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}$$


It sounds like you want inpainting, i.e., filling in missing parts of an image as best as possible (consistent with a cost function). There's a tremendous amount of research that's been done on schemes for inpainting. I suggest you start by reading up on the standard schemes for this. See also Photoshop's content-aware fill feature, which is another name for a similar functionality.

So, rather than trying to invent your own method for inpainting, you might want to try using an existing method for inpainting (for filling in the non-anchor pixel values). Since the inpainting problem has been studied so extensively in the literature, existing methods might perform better than anything you can invent in a few weeks.

• Yes, it can also be used for Inpainting. Yet I want this method. Do you know how to solve it? Should I post it on StackOverflow? – Royi Apr 15 '14 at 18:45
• @Drazick, what I'm suggesting is that you might be able to use an existing inpainting method to fill in the non-anchor pixel values (rather than trying to invent your own method). – D.W. Apr 16 '14 at 0:38
• Have a look above, I managed to proceed. Do you think you can assist? – Royi Apr 21 '14 at 22:37