{ "cells": [ { "cell_type": "markdown", "id": "77eb0b9a", "metadata": {}, "source": [ "# Lesson 9: Image Warping, Inverse Mapping, and Bilinear Interpolation\n", "\n", "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*." ] }, { "cell_type": "code", "execution_count": null, "id": "7ce89699", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import cv2\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "id": "bc4e3ed9", "metadata": {}, "outputs": [], "source": [ "def make_test_image(size=60):\n", " img = np.zeros((size, size, 3), dtype=np.uint8)\n", " cv2.rectangle(img, (5, 5), (size - 5, size - 5), (60, 90, 160), -1) # muted blue\n", " tri = np.array([[30, 12], [10, 32], [50, 32]], dtype=np.int32)\n", " cv2.fillPoly(img, [tri], (225, 195, 60)) # muted yellow\n", " return img\n", "\n", "img = make_test_image()\n", "plt.imshow(img)\n", "plt.title('Test image')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ea9f5810", "metadata": {}, "source": [ "## The problem with forward mapping\n", "\n", "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**.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "3619afe5", "metadata": {}, "outputs": [], "source": [ "h, w = img.shape[:2]\n", "M = cv2.getRotationMatrix2D((w / 2, h / 2), angle=25, scale=1.6)\n", "\n", "forward = np.zeros_like(img)\n", "ys, xs = np.mgrid[0:h, 0:w]\n", "src_pts = np.stack([xs.ravel(), ys.ravel(), np.ones(xs.size)])\n", "dst_pts = M @ src_pts\n", "\n", "dxi = np.round(dst_pts[0]).astype(int)\n", "dyi = np.round(dst_pts[1]).astype(int)\n", "valid = (dxi >= 0) & (dxi < w) & (dyi >= 0) & (dyi < h)\n", "forward[dyi[valid], dxi[valid]] = img[ys.ravel()[valid], xs.ravel()[valid]]\n", "\n", "filled = np.count_nonzero(forward.sum(axis=2))\n", "print(f'{filled} / {h*w} destination pixels received a value ({100*filled/(h*w):.0f}%)')\n", "\n", "plt.imshow(forward)\n", "plt.title('Forward mapping: visible holes (black speckle)')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "1e80d62f", "metadata": {}, "source": [ "## The fix: inverse mapping\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "b36426f3", "metadata": {}, "outputs": [], "source": [ "M_inv = cv2.invertAffineTransform(M)\n", "\n", "dst_grid = np.stack([xs.ravel(), ys.ravel(), np.ones(xs.size)])\n", "src_coords = M_inv @ dst_grid\n", "src_x = src_coords[0].reshape(h, w)\n", "src_y = src_coords[1].reshape(h, w)\n", "\n", "print('example fractional source coordinate for destination pixel (10, 15):')\n", "print(f' ({src_x[15, 10]:.3f}, {src_y[15, 10]:.3f})')" ] }, { "cell_type": "markdown", "id": "f87c9897", "metadata": {}, "source": [ "### Nearest-neighbor interpolation\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "fcab0334", "metadata": {}, "outputs": [], "source": [ "def nearest_sample(image, xf, yf):\n", " h, w = image.shape[:2]\n", " xi = np.round(xf).astype(int)\n", " yi = np.round(yf).astype(int)\n", " valid = (xi >= 0) & (xi < w) & (yi >= 0) & (yi < h)\n", " out = np.zeros(xf.shape + (image.shape[2],), dtype=np.uint8)\n", " out[valid] = image[yi[valid], xi[valid]]\n", " return out\n", "\n", "nearest = nearest_sample(img, src_x, src_y)\n", "\n", "plt.imshow(nearest)\n", "plt.title('Inverse mapping + nearest neighbor: no holes, but blocky')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "5d317d6f", "metadata": {}, "source": [ "### Bilinear interpolation\n", "\n", "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:\n", "\n", "$$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)$$\n", "\n", "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 *bi*linear." ] }, { "cell_type": "code", "execution_count": null, "id": "34929721", "metadata": {}, "outputs": [], "source": [ "def bilinear_sample(image, xf, yf):\n", " h, w = image.shape[:2]\n", " x0 = np.floor(xf).astype(int)\n", " y0 = np.floor(yf).astype(int)\n", " x1, y1 = x0 + 1, y0 + 1\n", " dx = (xf - x0)[..., None]\n", " dy = (yf - y0)[..., None]\n", "\n", " valid = (x0 >= 0) & (x1 < w) & (y0 >= 0) & (y1 < h)\n", " x0c, x1c = np.clip(x0, 0, w - 1), np.clip(x1, 0, w - 1)\n", " y0c, y1c = np.clip(y0, 0, h - 1), np.clip(y1, 0, h - 1)\n", "\n", " Ia = image[y0c, x0c].astype(np.float64)\n", " Ib = image[y0c, x1c].astype(np.float64)\n", " Ic = image[y1c, x0c].astype(np.float64)\n", " Id = image[y1c, x1c].astype(np.float64)\n", "\n", " out = (1 - dx) * (1 - dy) * Ia + dx * (1 - dy) * Ib + (1 - dx) * dy * Ic + dx * dy * Id\n", " out[~valid] = 0\n", " return np.clip(out, 0, 255).astype(np.uint8)\n", "\n", "bilinear = bilinear_sample(img, src_x, src_y)\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))\n", "axes[0].imshow(nearest)\n", "axes[0].set_title('Nearest neighbor')\n", "axes[0].axis('off')\n", "axes[1].imshow(bilinear)\n", "axes[1].set_title('Bilinear')\n", "axes[1].axis('off')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "bfb1cc99", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "af43a725", "metadata": {}, "source": [ "### Sanity check against OpenCV\n", "\n", "`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)." ] }, { "cell_type": "code", "execution_count": null, "id": "7e0f9c93", "metadata": {}, "outputs": [], "source": [ "cv_result = cv2.warpAffine(img, M, (w, h), flags=cv2.INTER_LINEAR)\n", "\n", "diff = np.abs(bilinear.astype(int) - cv_result.astype(int))\n", "print(f'max pixel difference = {diff.max()} (out of 255)')\n", "print(f'mean pixel difference = {diff.mean():.3f}')\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(9, 3.5))\n", "for ax, im, title in zip(axes, [bilinear, cv_result, diff.astype(np.uint8) * 20],\n", " ['Our bilinear', 'cv2.warpAffine', 'Difference (x20)']):\n", " ax.imshow(im)\n", " ax.set_title(title, fontsize=9)\n", " ax.axis('off')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "7fc01658", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "a63d87f0", "metadata": {}, "source": [ "### Exercise\n", "\n", "1. 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.\n", "2. 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?\n", "3. 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?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }