Interpolation Optimization Problem

Question

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.

Answer 1

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} $$

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 the methods to solve this problem by blocks to avoid memory limits?

Answer 2

$\newcommand{\be}{\mathbf{e}}\newcommand{\bL}{\mathbf{L}}\newcommand{\bC}{\mathbf{C}}\newcommand{\bt}{\mathbf{t}}\newcommand{\bl}{\mathbf{λ}}$Given Roy's formulation, we can solve this exactly without the need to invoke a QP solver.

To see why, we can add the constraint into the objective function using a Lagrange multiplier. Formally, let $$f(\be, \bl) = \be' \bL \be + \bl'(\bC\be - \bt)$$ be the cost function for a given solution $\be, \bl$ where $\mathbf{x}'$ indicates the transpose of $\mathbf{x}$. Dispite this being a quadratic function of $\be$ and linear function of $\bl$ there exists a unique solution that can be found by solving a system of linear equations.

To see why, let's first compute the gradients. The gradients of our unknowns are given by$$\begin{align*}\frac{\partial f}{\partial \be} &= 2\bL\be + \bC'\bl\\ \frac{\partial f}{\partial \bl} &= \bC\be - \bt.\end{align*}$$ Now to optimize for each parameter we see this is another system of linear equations, $$\begin{bmatrix}2\bL & \bC' \\ \bC & \mathbf{0} \end{bmatrix}\begin{bmatrix}\be \\ \bl \end{bmatrix} = \begin{bmatrix}\mathbf{0} \\ \bt \end{bmatrix}$$

The unique solution (provided that $\bL$ and $\bC\bL^{-1}\bC'$ are non-singular) for $\be$ and $\bl$ are determined by $$\begin{align*}\hat{\be} &= -\frac{1}{2} \bL^{-1}\bC'\bl\\ \hat{\bl} &= -2(\bC \bL^{-1} \bC')^{-1}\bt.\end{align*}$$

Solving this system in practice may still be an issue, however modern linear algebra libraries, such as BLAS or MKL, have built-in parallelization that can drastically speed things up.

Answer 3

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.

See also https://stackoverflow.com/q/8645515/781723 and https://stackoverflow.com/q/11338164/781723 and https://stackoverflow.com/q/2530449/781723

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.