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.

  • $\begingroup$ Have you tried using a QP solver like CPLEX? It seems you could restate this as a min-cost flow problem. $\endgroup$ Apr 15, 2014 at 2:16
  • $\begingroup$ Could you elaborate on that? I am not familiar with it. $\endgroup$
    – Royi
    Apr 15, 2014 at 2:18
  • $\begingroup$ 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. $\endgroup$ Apr 15, 2014 at 2:43
  • $\begingroup$ What do you mean by "expected results"; can you show an example? Are you sure the implementation and not the formalization is the issue? $\endgroup$
    – Raphael
    Apr 15, 2014 at 10:22
  • $\begingroup$ 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. $\endgroup$
    – Royi
    Apr 15, 2014 at 11:08

3 Answers 3


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?
  • $\begingroup$ Feel free to edit and add information to my answer as it is not complete (How should the sparse matrix be built efficiently? How can we calculate specifically $ {L}_{ij} $?). $\endgroup$
    – Royi
    Jul 18, 2014 at 5:20

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

  • $\begingroup$ I think the matrix $(I - A)$ is SDD (if not Laplacian) and very fast (in theory) solvers may be utilized to exploit the combinatorial structure of the underlying graph (e.g., sites.google.com/a/yale.edu/laplacian). $\endgroup$ May 3, 2017 at 18:01

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.

  • $\begingroup$ 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? $\endgroup$
    – Royi
    Apr 15, 2014 at 18:45
  • $\begingroup$ @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). $\endgroup$
    – D.W.
    Apr 16, 2014 at 0:38
  • $\begingroup$ Have a look above, I managed to proceed. Do you think you can assist? $\endgroup$
    – Royi
    Apr 21, 2014 at 22:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.