Many astronomers, particularly in imaging-focused research areas like my own, have used image alignment code written by others without struggling through writing their own.

Sometimes, you really should just use the library code so you can publish in a timely fashion. Really, though, you should understand exactly what's going on. Otherwise, you may find yourself with pointed questions from your coauthors (or your referees) that you're unable to answer!

Fast image interpolation is a cornerstone of computer graphics. We take this for granted today, but being able to zoom, scale, translate, and rotate images requires doing operations on thousands of pixels at sixty frames-per-second for applications like image editors or animation tools.

Since it is so fundamental, this topic has already been explained many ways by many people. Now I'm adding my own version to the mix.

The preliminaries¶

First, we need the usual Python tools for numerical computing and plotting:

Make the plots a little more science-y:

The first order of business is to get an image to interpolate:

That makes an image like this.

The "F" shape will make it clear when we flip the parity of the image, as well as showing which direction we're rotating when we rotate.

The simplest transformation¶

The simplest transformation is a full-pixel (integer) translation. Let's make an empty image...

...and fill it with the above pattern, shifted ten pixels to the right.

(We're cheating a little and just using NumPy indexing behavior to handle out-of-bounds coordinates—i.e. a negative src_x. Moving ten pixels left wouldn't work.)

There are a few things that should stand out to you about the above code.

Thing 1: loops¶

The first is that we have loops: in fact, a nested loop, which implies $O(n^2)$ scaling in the dimension of the image. There's no way around this; we have to visit every pixel to set its value. However, Python users know loops are slow-ish.

Let's compare a simple nested loop that sums all the pixels in the src image. First, the Python implementation:

Not so bad, unless you do it a few thousand times. But it could be much faster:

NumPy is still looping over all the elements of the array, but it does so in optimized C code instead of Python. That's enough to get a three orders of magnitude speedup.

So, be warned: the interpolation code developed in this post is optimized for clarity, not speed. (Maybe in another blog post!)

Thing 2: loop over dest, not src¶

It was unintuitive (at least, to me) that the interpolation function to transform from some image src to some image dest involves looping over the pixels in dest. But, of course, the destination's pixel coordinates are fixed. We can't write a value into the 2.75'th pixel on an axis, just the 2nd or the 3rd pixel. So, we start with the destination pixel coordinates dest_x and dest_y and transform the coordinate to find the relevant part of the source image (that's the src_y, src_x = dest_y, dest_x - dx in the example above).

Coordinate transformations¶

The typical representation of coordinate transformations for 2D graphics is as an affine transformation matrix. For the coordinates dest_x, dest_y we define a coordinate vector $\vec{v} = \begin{bmatrix}x \\ y \\ 1\end{bmatrix}$. The extra $1$ entry makes this an "augmented" vector, meaning the transform matrix becomes a $3 \times 3$ matrix instead of $2 \times 2$.

Without this extra entry, the transformation matrix could not represent a translation of the original coordinates. To get the 2D coordinates after transforming, simply drop the last entry of the vector.

As a quick refresher, here are the affine transformation matrices. (See Wikipedia for more details, pictures, etc.)

Translation¶

$$ \begin{bmatrix}x + \Delta x \\ y + \Delta y \\ 1\end{bmatrix} = \begin{bmatrix} 1 & 0 & \Delta x\\ 0 & 1 & \Delta y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}x \\ y \\ 1\end{bmatrix} $$

Scaling¶

$$ \begin{bmatrix}s_x x\\ s_y y\\ 1\end{bmatrix} = \begin{bmatrix} s_x & 0 & 0\\ 0 & s_y & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}x \\ y \\ 1\end{bmatrix} $$

Note that a parity flip (i.e. reversing the X or Y direction) is just a special case of scaling with a scale factor < 0.

Rotation¶

$$ \begin{bmatrix}x'\\ y'\\ 1\end{bmatrix} = \begin{bmatrix} \cos \theta & -\sin \theta & 0\\ \sin \theta & \cos \theta & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}x \\ y \\ 1\end{bmatrix} $$

Shearing¶

$$ \begin{bmatrix}x'\\ y'\\ 1\end{bmatrix} = \begin{bmatrix} 1 & c_x & 0\\ c_y & 1 & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}x \\ y \\ 1\end{bmatrix} $$

Who's getting transformed, anyway?¶

To find the relevant part of src given dest_x, dest_y, we transform from the destination image's pixel coordinates back to source image pixel coordinates. That means we have to think carefully about what's moving.

Suppose you have a feature at $(1, 1)$ in the source image and you want to place it at $(2, 3)$ in the destination image (i.e. a displacement of $+1$ in $x$ and $+2$ in $y$). That means when you transform $(2, 3)$ you should get $(1, 1)$, which means the $\Delta x$ and $\Delta y$ in the example above should be negative:

$$ \begin{bmatrix} 1 & 0 & -1\\ 0 & 1 & -2\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}2 \\ 3 \\ 1\end{bmatrix} = \begin{bmatrix}2 - 1 \\ 3 - 2 \\ 1\end{bmatrix} =\begin{bmatrix}1 \\1\\ 1\end{bmatrix} $$

So, when constructing a transformation, you must always be clear about whether you're expressing the image's motion or the coordinate transformation "back" to the source image pixels.

Nearest-neighbor interpolation¶

Technically, interpolation just means estimating values for new points at which we have no data. Using the nearest point where we do have data is not very accurate, but for images with large regions of uniform values it does okay. It is also very easy to implement.

Taking our image transformation loop from above, we change the way src_x and src_y are computed:

So far, so good. But we're moving by whole pixels still, which is easy mode. What if we wanted to rotate by $\pi / 6$?

That's because we're no longer using only whole integers, so this is where the rounding to nearest neighboring pixels comes in. Let's just round().

Still not good enough. What's happening is obvious, but what to do about it is subtle. Not only do we not have values at exactly the positions in the destination pixel grid, but we also are trying to interpolate values from positions that are more than a pixel outside the region covered by the source image.

In other words, there are pixels we can interpolate, and pixels we have to extrapolate.

To keep things simple and make sure we have a valid src_x and src_y for every destination pixel, we can just clamp the source coordinates to the range [0, n_pix - 1].

Well, that got us past the exceptions. But there's the matter of whether you want the image to rotate by $\pi / 6$, and which direction should be positive. Conventionally, rotation counterclockwise from the $+x$ towards the $+y$ axis is positive. So, again, we negate the value in the transformation matrix:

Better! Let's wrap that into a function.

Looks the same, which we expected. But, if you expected the F to rotate nicely about the center, you're probably disappointed. The explanation is hinted at by the ticks on the coordinate axes: the origin is pixel (0, 0), so we're pivoting about that.

For my applications, it's more useful to have a rotation about the center. Let's construct one. First, some helper functions so we don't have to write out the transform matrices by hand.

Now we can perform a rotation about the center as a combination of three steps: move center to origin, perform rotation about origin, move origin to image center.

Note that we want to support destination images that aren't the same dimensions as the source image, so we need the shapes of both the source and destination images. (The translation amounts are different for the centers of the source and the destination images.)

Step 1: move the image's center to the origin¶

Equivalently, transform the pixel coordinates of the destination such that (0, 0) maps to the source image's center.

Step 2: apply rotation¶

Here we see the power of composable transformations through matrix products. By multiplying the affine transformation matrices together we can chain together image operations. However, note that the order matters, and we're transforming coordinates and not images. That means for our destination pixels to be a rotated view of the source, we de- rotate first, then translate.

Step 3: move the origin back to the center of the image¶

Since we have a dest_shape that's double the size of the source image, we translate to its center rather than simply undoing the original translation.

Then, we can put this all together in a helper function.

The same applies to scaling, so let's make a function that does both

Using the destination shape argument we can upscale images to a higher resolution:

Bilinear: more sophisticated interpolation¶

I said above that nearest-neighbor interpolation works best when there are large regions of the same value, which is certainly true of the F test image. Let's see a less-cooperative example. Here I have plotted $f(r) = \sin(5 \pi r)$.

The underlying function is clearly radially symmetric, so rotating it shouldn't change it--except that we've discretized it on a grid of pixels.

For fun, let's rotate it through $2\pi$ in angle and see if it differs from the input.

Well, that's unsastisfactory. We need to level up our interpolation function.

The simplest nontrivial interpolation function is bilinear, which takes the values at the four nearest pixels and interpolates three times.

Let us suppose (without loss of generality) there are only four pixels in the image, and we're trying to interpolate a value in between them. Then for the hypothetical continuous function $f(x, y)$ we only have known values for $f(0, 0)$, $f(0, 1)$, $f(1, 0)$, and $f(1, 1)$. We are looking for an approximate value for $f(\Delta x, \Delta y)$, where $0 \le \Delta x \le 1$ and $0 \le \Delta y \le 1$.

In the figure below, the known values at the source image pixel locations are shown by the red circular markers. Using the pair of points $f(0, 0)$ and $f(0, 1)$ to define a line, we can interpolate an approximate value for $f(\Delta x, 0)$ by assuming it lies on this line.

Linearly interpolating gives $f(\Delta x, 0) \approx f(0, 0) (1 - \Delta x) + f(1, 0)(\Delta x)$, the value at the lower square marker.

Repeating the procedure with $f(1, 0)$ and $f(1, 1)$ gives us an approximation for $f(\Delta x, 1) \approx f(1, 0) (1 - \Delta x) + f(1, 1)(\Delta x)$, shown by the upper square marker.

Finally, we interpolate between $f(\Delta x, 0)$ and $f(\Delta x, \Delta y)$ to get the interpolated value for our destination pixel: \begin{align*} f(\Delta x, \Delta y) \approx&\ f(\Delta x, 0) (1 - \Delta y) + f(\Delta x, 1)(\Delta y)\\ =&\ (f(0, 0) (1 - \Delta x) + f(1, 0)(\Delta x)) (1 - \Delta y) + (f(1, 0) (1 - \Delta x) + f(1, 1)(\Delta x))(\Delta y)\\ =&\ f(0, 0) (1 - \Delta x) (1 - \Delta y)\\ &\ +\ f(1, 0) (\Delta x) (1 - \Delta y)\\ &\ +\ f(1, 0) (1 - \Delta x) (\Delta y)\\ &\ +\ f(1, 1)(\Delta x)(\Delta y) \end{align*}

We only work this out for the case where we're in the first pixel in the grid because we can always choose our coordinates such that we are there. Let's make a function to do this calculation, cutting out the $2 \times 2$ pixels we care about from the source image.

The comments point out the places where we handle out-of-bounds coordinates by returning not-a-number (a.k.a. NaN or np.nan).

Updating the simple_image_xform function to use bilinear_patch to find the destination pixel values, we have bilinear_image_xform:

Two things (again) to notice here:

Improved residuals¶

We can tell we're doing better because the residuals (rightmost panel) are smaller on average than in the nearest-neighbor example above. So, even though the image looks different up close, we're getting closer to the ideal (i.e. rotational symmetry).

Missing pixels¶

From the very first frame, there is one row and one column that are white (i.e. missing). The destination pixel no longer depends on just one source pixel; we need a little domain of four pixels.

On frame zero, we're rotating by 0º, which shouldn't change anything. However, for an $N \times N$ image, the corner pixel at $(N - 1, N - 1)$ will depend on values in pixels $(N - 1, N)$, $(N, N - 1)$, and $(N, N)$. Of course, those values are unknown (or unknowable), so the result of the interpolation is undefined.

There are various reasonable workarounds, like checking if the transformation maps to an integer coordinate in the source image and skipping the whole interpolation step. Another option is to replace the missing values in the bilinear_patch calculation with the nearest known values. Both are justifiable, but you should know which one you're doing so that you don't get confused by interpolation artifacts.

Bicubic interpolation or cubic convolution: best speed/accuracy tradeoff¶

Bilinear interpolation assumes the underlying continuous function sampled on our pixel grid increases or decreases only linearly in the space between pixels. Of course, we know plenty of functions without constant first derivatives. Can we capture some of this higher order behavior? The answer is yes, at the cost of some complexity and additional arithmetic operations.

With a carefully chosen kernel function, a $4 \times 4$ domain of pixels can contribute to a cubic curve that connects two of the pixels. (This is not quite spline interpolation, but a good approximation that doesn't require derivative information for the underlying function.)

To implement bicubic interpolation let's start with an example image. The mini_img below is a $4 \times 4$ image covering the region $-1 \le x \le 2$ and $-1 \le y \le 2$.

(Nearest-neighbor uses one pixel, bilinear uses $2 \times 2$ pixels, and bicubic uses $4 \times 4$. Why doesn't anyone use $3 \times 3$, for biquadratic? Well... that's kind of complicated.)

If we use the 1D cubic convolution kernel to interpolate values between columns, we get four continuous curves representing $f(x, -1)$, $f(x, 0)$, $f(x, 1)$, $f(x, 2)$, one for each row.

As above with bilinear interpolation, we are assuming (without loss of generality) that we are interpolating coordinates in the domain $0 < \Delta x < 1$ and $0 < \Delta y < 1$. Evaluating these approximations at $x = \Delta x$ for each gives us four new values, $f(\Delta x, -1)$, $f(\Delta y, 0)$, $f(\Delta y, 1)$, $f(\Delta y, 2)$.

Plugging these into the same 1D interpolation function allows us to approximate $f(\Delta x, \Delta y)$.

This (somewhat complicated) figure shows what's going on.

First, each row is interpolated at $\Delta x$, shown as smaller circular points between 0 and 1. The curve connecting $y = 0$ to $y = 1$ is then evaluated at $y = \Delta y$ to get an approximate value for $f(\Delta x, \Delta y)$.

What's the magic kernel function, then? It was derived by one Robert G. Keys in a 1981 publication ("Cubic convolution interpolation for digital image processing", IEEE Transactions on Acoustics, Speech, and Signal Processing). It's not even that hard to follow if you have some knowledge of Taylor series and algebra.

(Unfortunately, as of this writing IEEE is still trying to charge $33.00 for access to a paper that's more than four decades old. So you'll have to find it on your own.)

The interpolation kernel is

$$ u(s) = \begin{cases} \frac{3}{2} |s|^3 - \frac{5}{2} |s|^2 + 1 & 0 < |s| < 1\\ -\frac{1}{2} |s|^3 + \frac{5}{2} |s|^2 - 4 |s| + 2 & 1 < |s| < 2\\ 0 & s < |s| \end{cases} $$

which, with some algebraic manipulation given in the paper, can produce the interpolation function $g(x) \approx f(x)$ which can be evaluated to approximate $f(x)$ between samples. The interpolation function is

\begin{align*} g(x) =&\ c_{k - 1} (- s^3 + 2 s^2 - s)/2 + c_k (3s^3 - 5s^2 + 2)/2\\ &+ c_{k + 1} (-3 s^3 + 4s^2 + s)/2 + c_{k + 2}(s^3 - s^2)/2 \end{align*}

where $s = (x - x_k) / h$, $h$ is the sampling interval, and $c_k = f(x_k)$ for $k = 0, 1, 2, ..., N$ are the known sample values at the coordinates $x_k$.

We've set up the coordinates above to run from -1 to 2, and we're using pixel coordinates so $h = 1$ and can be neglected. Therefore $s = \Delta x$ and $c_{k - 1} = f(-1)$, $c_{k} = f(0)$, and so on.

(Note: Here we're making a slight abuse of notation, since we're interpolating in 1D: $f(x)$ is really $f(x, y)$ for some sampled y coordinate.)

The simplified expression for our chosen coordinates becomes

\begin{align*} g(x) = &f(-1) \frac{-(\Delta x)^3 + 2(\Delta x)^2 - (\Delta x)}{2}\\ &+ f(0) \frac{3(\Delta x)^3 - 5(\Delta x)^2 + 2}{2}\\ &+ f(1) \frac{-3(\Delta x)^3 + 4(\Delta x)^2 + \Delta x}{2}\\ &+ f(2) \frac{(\Delta x)^3 - (\Delta x)^2}{2} \end{align*}

We then factor out the 1/2 and group terms by powers of $\Delta x$:

\begin{align*} g(x) = \frac{1}{2} [ &\left[-f(-1) + 3f(0) - 3f(1) + f(2)\right] (\Delta x)^3 \\ &+ \left[ 2f(-1) - 5f(0) + 4f(1) - f(2)\right] (\Delta x)^2 \\ &+ \left[ -f(-1) + f(1) \right] (\Delta x)\\ &+ 2f(0) ] \end{align*}

which we can translate to Python code:

Now let's interpolate $f(\Delta x, \Delta y)$ in the image above for some values of $\Delta x$ and $\Delta y$.

The interpolation curve from cubic1d should pass through $f(0, y)$ and $f(1, y)$ when we evaluate it on the interval $0 < x < 1$.

This shows the interpolated value for $f(\Delta x, 0)$ as an orange dot on the curve. (Why not draw a curve through all three pairs? Well, you need a domain of $4 \times 4$ pixels, so we'd need one more pixel on the "low" side and two more on the "high" side.)

To make this 2D we repeat for every row of the image:

Above, we've plotted the values from the different rows offset slightly for clarity.

Taking the interpolated values for $f(\Delta x, y)$ for $y = -1, 0, 1, 2$ (shown by the point markers on the interpolated curves above) as a new 1D sequence, we interpolate again:

This interpolation along a second direction is the "bi" in bicubic. We can write it up as a function that assumes a 4 x 4 pixel subimage to operate on.

The interpolated $f(\Delta x, \Delta y)$ should agree, and it does:

Now, we just have to do this for every pixel in the destination image. Easy!

Since it requires a $4 \times 4$ pixel domain to interpolate from, I've restructured the code to construct a cutout image with a fill value for missing pixels (i.e. on the edges).

Let's see the effect of different fill values on the resulting image. Here I'm rotating by 45º with zeros and with np.nans for out-of-bounds points.

Arguably, the interpolated values at the edges on the left are just nonsense and there for aesthetic reasons. The missing value treatment least likely to mislead you is to treat every computation that takes an undefined value as producing an undefined value, which is exactly what a fill value of np.nan does.

Let's see how it does with the radially symmetric signal.

Not bad! Returning to the "F" test image for a moment, let's investigte the bicubic interpolator's behavior around sharp edges. Here we translate the test pattern by a noninteger number of pixels (to make the sub-pixel behavior of the interpolator visible).

Let's zoom in on the base of the F and plot a slice through the image.

We put in ones and zeros in the test pattern, so why have we got values below zero and above one now? In the test pattern, the the edge of the "F" in the test image is a discontinuity; from one pixel to the next the value jumps from 0 to 1. Intuitively, we think that the edge should be "infinitely sharp", so the result of translation should still transition from zero to one in the space of a single pixel.

Why is this so different from our intuition? Is this "right"?

Oversampling, undersampling, and band limits¶

This is "ringing", basically the same thing as Gibbs phenomenon. It all boils down to what you're capturing in your pixels. If your image contains infinitely sharp edges, like our "F" example, you need an infinite number of pixels to express that in a pixelated / sampled format. So, no, the image above is not "right"; this is why graphic designers use a completely different image representation to work in. In practice, we have given the computer a number of samples, representing a signal up to some band-limit that corresponds to how finely sampled (i.e. sharp) it is.

Fortunately, if you're taking pictures of natural phenomena, there are no infinitely sharp edges. When considering the detection of exoplanet signals, we are dealing with a more constrained version of the image interpolation problem than, e.g., Photoshop does. In high-contrast imaging, our science images are always oversampled relative to the band-limit of the telescope system. That means that we can be confident that our interpolation and resampling is preserving information (at least, to a good approximation).

Consider a counter-example: an undersampled image of a star might look like this:

Here I've marked the star location at $(2, 2)$. It might really be at $(2.25, 2)$, but it's such a tiny feature it gets summed up in the pixel ranging from $1.5 < x < 2.5$ and $1.5 < y < 2.5$.

What if we translate the image 0.5 pixels right, using our interpolator from above?

If your model for the "star" in pixel $(2, 2)$ is that it's illuminating one pixel-sized square uniformly from its center, shifting the square by half a pixel puts the star on the border between pixels at $(2.5, 2)$, so it is quite reasonable to put half the light in the adjacent pixel. On the other hand, if your star is really at $(2.25, 2.25)$, the star shifts to $(2.75, 2.25)$ so more than half of the light should be there.

Given this transformed image, you would also be unable to tell if it is one star of brightness 1.0 on a pixel boundary, or two stars of brightness 0.5 that happen to be next to each other.

In an undersampled image, there's no way to know. In other words, there is no way to rotate or translate the undersampled image without introducing potentially major changes in it.

Though it's not an isolated point like a star, the hard edge of the "F" pattern is an example of one of these undersampled features.

Oversampled scene¶

The most extreme case for an undersampled scene is a single pixel, but what model should we use for an oversampled scene? You won't go far wrong with a 2D Gaussian; even extreme adaptive optics PSFs are nearly Gaussian in their core.

Let's construct a 16 x 16 pixel scene with a Gaussian at the center. To compare with the interpolated image we will also generate one where the Gaussian has been displaced by dx = 2.5 pixels. Here we're generating two "truth" images based on the true analytical expression for $f(x, y)$ (implemented as gauss2d).

Now we attempt to perform the translation by interpolation, using the "before" image.

That is close enough to the "truth" image that it's probably easier to make sense of a plot showing their difference.

Here you can see that the bicubic interpolation is not perfect. It is not possible in general to perfectly describe a Gaussian curve with a 3rd order polynomial, so you can see deviations where the cubic approximation is insufficient.

The right way to interpolate¶

Hopefully, the above has demonstrated that there is no one "right way" to resample images. Every method comes with trade-offs in speed and accuracy, and there are good ones I didn't cover (because this post is long enough). For a true band-limited signal, for example, using a zero-padded Fourier transform upsamples without loss of information (except for floating-point roundoff error). In practice, features like hot pixels and cosmic rays will always be "undersampled", and masking bad pixels will introduce ringing at the mask edges, so I don't use that method very often. Maybe I'll write that one up later.

The implementations above are intended to illuminate the way the algorithm works, but are not optimized for speed. Using someone else's optimized implementation is preferable if you don't have the time to build your own, but you should be sure to understand the accuracy and tradeoffs made by the implementation you use. For examples of the possible pitfalls, see the Bicubic bake-off post.