Lesson 9: Image Warping, Inverse Mapping, and Bilinear Interpolation¶

Lesson 8 built transformation matrices; this lesson is about actually applying one to resample an image. The naive way to do this — push every source pixel to its transformed location — turns out to be broken. We'll see why, fix it with inverse mapping, and then deal with the fact that inverse mapping lands on fractional pixel coordinates that need to be interpolated.

In [1]:

import numpy as np
import cv2
import matplotlib.pyplot as plt

In [2]:

def make_test_image(size=60):
    img = np.zeros((size, size, 3), dtype=np.uint8)
    cv2.rectangle(img, (5, 5), (size - 5, size - 5), (60, 90, 160), -1)   # muted blue
    tri = np.array([[30, 12], [10, 32], [50, 32]], dtype=np.int32)
    cv2.fillPoly(img, [tri], (225, 195, 60))  # muted yellow
    return img

img = make_test_image()
plt.imshow(img)
plt.title('Test image')
plt.axis('off')
plt.show()

$No description has been provided for this image$

The problem with forward mapping¶

The obvious way to warp an image: for every source pixel $(x,y)$, compute its destination $(x',y')=M(x,y)$ and copy the pixel value there. This is called forward mapping.

It has a fundamental flaw: whereas the source pixels form a complete grid, their transformed locations generally do not. When a transform enlarges the image (or rotates it, which locally stretches the grid along the diagonal), the destination locations spread out and leave gaps — no source pixel happens to land exactly on some destination pixels.

In [3]:

h, w = img.shape[:2]
M = cv2.getRotationMatrix2D((w / 2, h / 2), angle=25, scale=1.6)

forward = np.zeros_like(img)
ys, xs = np.mgrid[0:h, 0:w]
src_pts = np.stack([xs.ravel(), ys.ravel(), np.ones(xs.size)])
dst_pts = M @ src_pts

dxi = np.round(dst_pts[0]).astype(int)
dyi = np.round(dst_pts[1]).astype(int)
valid = (dxi >= 0) & (dxi < w) & (dyi >= 0) & (dyi < h)
forward[dyi[valid], dxi[valid]] = img[ys.ravel()[valid], xs.ravel()[valid]]

filled = np.count_nonzero(forward.sum(axis=2))
print(f'{filled} / {h*w} destination pixels received a value ({100*filled/(h*w):.0f}%)')

plt.imshow(forward)
plt.title('Forward mapping: visible holes (black speckle)')
plt.axis('off')
plt.show()

1408 / 3600 destination pixels received a value (39%)

No description has been provided for this image

The fix: inverse mapping¶

Instead of asking "where does each source pixel go?", ask the opposite question for every destination pixel: "where in the source image did this come from?" That means applying the inverse transform $M^{-1}$ to each destination coordinate. Since we now iterate over a complete destination grid, every output pixel is guaranteed to get a value — no holes, by construction.

The catch: $M^{-1}(x',y')$ almost never lands exactly on an integer source coordinate. We need to interpolate a value from the surrounding source pixels.

In [4]:

M_inv = cv2.invertAffineTransform(M)

dst_grid = np.stack([xs.ravel(), ys.ravel(), np.ones(xs.size)])
src_coords = M_inv @ dst_grid
src_x = src_coords[0].reshape(h, w)
src_y = src_coords[1].reshape(h, w)

print('example fractional source coordinate for destination pixel (10, 15):')
print(f'  ({src_x[15, 10]:.3f}, {src_y[15, 10]:.3f})')

example fractional source coordinate for destination pixel (10, 15):
  (22.633, 16.221)

Nearest-neighbor interpolation¶

The simplest option: just round to the nearest integer source pixel. This has no holes, but it's blocky — many destination pixels round to the same source pixel, and the image looks pixelated wherever the transform enlarges the image.

In [5]:

def nearest_sample(image, xf, yf):
    h, w = image.shape[:2]
    xi = np.round(xf).astype(int)
    yi = np.round(yf).astype(int)
    valid = (xi >= 0) & (xi < w) & (yi >= 0) & (yi < h)
    out = np.zeros(xf.shape + (image.shape[2],), dtype=np.uint8)
    out[valid] = image[yi[valid], xi[valid]]
    return out

nearest = nearest_sample(img, src_x, src_y)

plt.imshow(nearest)
plt.title('Inverse mapping + nearest neighbor: no holes, but blocky')
plt.axis('off')
plt.show()

Bilinear interpolation¶

Instead of snapping to the single nearest pixel, bilinear interpolation blends the 4 neighboring pixels, weighted by how close it is to each one. Let $(x_0,y_0)$ be the integer pixel just below-left of $(x,y)$, and let $d_x = x-x_0$, $d_y=y-y_0$, where $0 \le d_x, dy < 1$ are fractional coordinates. Then the interpolated value is:

$$I(x,y) \approx (1-d_x)(1-d_y)\,I(x_0,y_0) + d_x(1-d_y)\,I(x_1,y_0) + (1-d_x)\,d_y\,I(x_0,y_1) + d_x\,d_y\,I(x_1,y_1)$$

where $x_1=x_0+1$, $y_1=y_0+1$. This is just two 1D linear interpolations (along $x$, then along $y$) composed — hence bilinear.

In [6]:

def bilinear_sample(image, xf, yf):
    h, w = image.shape[:2]
    x0 = np.floor(xf).astype(int)
    y0 = np.floor(yf).astype(int)
    x1, y1 = x0 + 1, y0 + 1
    dx = (xf - x0)[..., None]
    dy = (yf - y0)[..., None]

    valid = (x0 >= 0) & (x1 < w) & (y0 >= 0) & (y1 < h)
    x0c, x1c = np.clip(x0, 0, w - 1), np.clip(x1, 0, w - 1)
    y0c, y1c = np.clip(y0, 0, h - 1), np.clip(y1, 0, h - 1)

    Ia = image[y0c, x0c].astype(np.float64)
    Ib = image[y0c, x1c].astype(np.float64)
    Ic = image[y1c, x0c].astype(np.float64)
    Id = image[y1c, x1c].astype(np.float64)

    out = (1 - dx) * (1 - dy) * Ia + dx * (1 - dy) * Ib + (1 - dx) * dy * Ic + dx * dy * Id
    out[~valid] = 0
    return np.clip(out, 0, 255).astype(np.uint8)

bilinear = bilinear_sample(img, src_x, src_y)

fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
axes[0].imshow(nearest)
axes[0].set_title('Nearest neighbor')
axes[0].axis('off')
axes[1].imshow(bilinear)
axes[1].set_title('Bilinear')
axes[1].axis('off')
plt.tight_layout()
plt.show()

The bilinear result has smooth edges around the triangle and rectangle instead of jagged steps — the classic trade-off is a softer, slightly blurrier image in exchange for removing aliasing artifacts.

Sanity check against OpenCV¶

cv2.warpAffine does exactly this — inverse mapping plus interpolation — internally. We compare our from-scratch version against cv2.warpAffine(..., flags=cv2.INTER_LINEAR) applied with the forward matrix M (OpenCV inverts it internally by default).

In [7]:

cv_result = cv2.warpAffine(img, M, (w, h), flags=cv2.INTER_LINEAR)

diff = np.abs(bilinear.astype(int) - cv_result.astype(int))
print(f'max pixel difference  = {diff.max()}   (out of 255)')
print(f'mean pixel difference = {diff.mean():.3f}')

fig, axes = plt.subplots(1, 3, figsize=(9, 3.5))
for ax, im, title in zip(axes, [bilinear, cv_result, diff.astype(np.uint8) * 20],
                          ['Our bilinear', 'cv2.warpAffine', 'Difference (x20)']):
    ax.imshow(im)
    ax.set_title(title, fontsize=9)
    ax.axis('off')
plt.tight_layout()
plt.show()

max pixel difference  = 1   (out of 255)
mean pixel difference = 0.092

The two match almost exactly — the tiny remaining difference (at most 1 gray level) comes from OpenCV using fixed-point rounding internally for speed, rather than full floating-point arithmetic.

Exercise¶

Why doesn't the forward-mapping holes problem happen when the transform only shrinks the image? Try a scale of 0.5 instead of 1.6 in the forward-mapping example and see what fraction of destination pixels get filled.
Extend bilinear_sample to bicubic interpolation conceptually: instead of a 2x2 neighborhood, what neighborhood size would you need, and what property would you want the weights to satisfy at the sample points?
Time your bilinear_sample against cv2.warpAffine on a larger image (e.g. 800x800) using %timeit. By how much faster is the OpenCV version, and why?

← Lesson 8: Geometric Transformations

Lesson 10: Convolution →