This semester I am taking two courses that are related to image and computer vision and I started to better understand that some amazing applications require even basic mathematical tools that we learned in the first two years in college. Yet many open source project now enables us to do lots of image analysis without a deep understanding of how basic mathematics work underneath, It is still fruitful and inspiring for me to all of sudden connect all the dots that we once collected.
Here, I want to introduce one of the examples – Image Interpolation.
What is image interpolation and why that is important?
It is convenient for us to zoom an image(say transform a \(500\times500\) image to \(1000\times1000)\) thanks to the help of the digital device. But when an image is zoomed up, there are actually extra pixels inserted into the original one, and how to color those extra pixels is a critical task called image interpolation.
A naive method
A reasonable image is always locally continuous. That’s saying that we are expecting neighboring pixels have similar colors. Thus, a naive way to fill one extra pixel is by copying the color from the nearest pixel. Noticing that there is still one more problem undefined – how to find the nearest pixel?
Let’s assume we have two pixel coordinates (see Fig1, one is \(s\) for the original image and another is \(S\) for the zoomed image. And the image is scaled by \(k\) times.
\(\forall (u, v) \in S\), the corresponding pixel in s is \((\frac{u}{k}, \frac{v}{k})\) and four candidate neighboring pixels in s to are: \((i, j), (i, j+1), (i+1, j), (i+1,j+1)\), where \(i = floor(\frac{u}{k})\) and \(j = floor(\frac{v}{k})\). And Define \(\Delta = (\Delta_i, \Delta_j) = (\frac{u}{k} – i, \frac{v}{k} – j)\)
Finally, the distance we are discussing here is the euclidean distance. Now we could find the nearest pixel to POI.
An improvement
Since we know how close each candidate pixel and my POI are, why not use a weighted average to gain a better estimation?
If possible, then how to select weights wisely? Intuitively, we can use the area as a set of weights (which requires fewer computations than euclidean distance and all the weights (total area) magically sum up to 1!).
Bilinear Image Interpolation
In another word, the technique we get from the improvement above is:
$$f(i+\Delta_i, j+\Delta_j) = \sum_{i, j} f(i,j)S_{i,j}$$
And it also has a fancy name – bilinear image interpolation.
Mathematical Foundation Behind
It idea of the weighted average is natural and intuitive, yet we can find a solid mathematical foundation from numerical analysis – Interpolation.
Considering one dimension, we have multiple ways to find the interpolation function giving a set of points. For instance, if we are given two points \((x_i, f_i)\), \((x_{i+1}, f_{i+1})\) then using linear polynomial we could estimate \(f(x), \forall x in [i, i+1]\). The Lagrange interpolation form can be written as:
$$f(x) = f_{i}\frac{x-x_{i+1}}{(x_{i} – x_{i+1})} + f_{i+1}\frac{x-x_{i}}{(x_{i+1} – x_{i})}$$
In two dimensions, we could use the same trick three times to estimate our POI (See Fig3).
- estimate \( f_{i, j+\Delta_{j}} \) based on \(f_{i, j} , f_{i,j+1}\)
- estimate \(f_{i+1, j+\Delta_j}\) based on \(f_{i+1, j} , f_{i+1,j+1}\)
- estimate \(f_{i+\Delta_i, j+\Delta_j}\) based on \(f_{i, j+\Delta_j} , f_{i+1,j+\Delta_j}\)
Consider using linear algebra form to represent the Lagrange interpolation process steps above:
$$f_{i, j+\Delta_j} =\begin{bmatrix}f_{i, j} & f_{i,j+1}\end{bmatrix}\begin{bmatrix}1-\Delta_j \\ \Delta_j\end{bmatrix}$$
$$f_{i+1, j+\Delta_j} =\begin{bmatrix}f_{i+1, j} & f_{i+1,j+1}\end{bmatrix}\begin{bmatrix}1-\Delta_j \\ \Delta_j\end{bmatrix}$$
$$f_{i+\Delta_i, j+\Delta_j} =\begin{bmatrix}1-\Delta_i & \Delta_i\end{bmatrix}\begin{bmatrix}f_{i+1, j+\Delta_j} \\ f_{i,j+\Delta_j}\end{bmatrix}$$
To sum up, the pixel value \(f_{i+\Delta_i, j+\Delta_j}\) can be written as:
\[f_{i+\Delta_i, j+\Delta_j} = \begin{bmatrix}1-\Delta_i & \Delta_i\end{bmatrix}\begin{bmatrix}f_{i,j} & f_{i,j+1} \\ f_{i+1, j} & f_{i+1, j+1} \end{bmatrix}\begin{bmatrix}1-\Delta_j \\ \Delta_j\end{bmatrix}\]
Noticing how amazing we could arrange those four pixels’ intensity in the middle matrix(say \(F\)) as if they are in the \(s\) coordinate system(See Fig4).
Now it is relatively easy to see that:
\[f_{i+\Delta_i, j+\Delta_j}=(1-\Delta_i)(1-\Delta_j)\cdot f_{i, j}+\Delta_{i}(1-\Delta_{j})\cdot f_{i, j+1}+(1-\Delta_{i})\Delta_{j}\cdot f_{i+1,j}+\Delta_{i}\Delta_{j}\cdot f_{i+1, j+1}\\=\sum_{i,j}S_{i,j}f(i,j)\]
which proves our intuitive idea at the beginning.
Bicubic Image Interpolation
Bicubic image interpolation is trying to consider more neighboring pixels \((4\times 4)\)when estimating our POI. Noticing the matrix form in the bilinear interpolation section above, without loss of generality, the universal model we are adapting is actually:
$$f_{i,j}=W(\Delta_j)^{T}FW(\Delta_i)$$
In one dimension, considering the Lagrange Interpolation function based on cubic polynomial given 4 points:
$$f(x)=\sum_{k=i-1}^{i+2}f_{k}\frac{\Pi_{l=i-1, l\neq k}^{i+2}(x-x_{l})}{\Pi_{l=i-1, l\neq k}^{i+2}(x_{k}-x_{l})}$$
Taking \(x=i+\Delta_i, x_{k} = k\) and simplifying the right-hand part of the equation above, we can see that the weight vector \(W(\Delta_i)\) should be:
\[\begin{bmatrix}-\frac{1}{6}\Delta_i(\Delta_i-1)(\Delta_i-2) \\ \frac{1}{2}(\Delta_i+1)(\Delta_i-1)(\Delta_i-2)\\-\frac{1}{2}(\Delta_i+1)\Delta_i(\Delta_i-2) \\ \frac{1}{6}(\Delta_i+1)\Delta_i(\Delta_i-1)\end{bmatrix}\]
Now, we have a clear and neat equation for estimating POI using the bicubic method.
Summary
This article shows a new way to understand image interpolation and mathematics from numerical analysis behind the scene. A unified matrix form of pixel interpolation is presented and also raised a feasible algorithm using the Lagrange interpolation formula to obtain weights for neighboring pixels.