{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3dcf84ac",
   "metadata": {},
   "source": [
    "# Lesson 10: Convolution\n",
    "\n",
    "Convolution is the core operation behind blurring, sharpening, and edge detection — and it's also the operation at the heart of a convolutional neural network's convolution layers (Lesson 33). This lesson builds it from scratch, clarifies a common point of confusion (convolution vs. correlation), and shows a few classic filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eb1fca1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d292b213",
   "metadata": {},
   "source": [
    "## 1D convolution\n",
    "\n",
    "Before tackling 2D images, let's first start with the simpler problem of convolving two 1D arrays. For a signal $I$ and a kernel $K$, the discrete convolution is\n",
    "\n",
    "$$(I * K)(x) = \\sum_{i} K(i)\\, I(x - i)$$\n",
    "\n",
    "Concretely: **flip** the kernel end-to-end, **slide** it along the signal one position at a time, and at each position take the dot product between the flipped kernel and the overlapping chunk of signal. Now let's try it with a signal that ramps up and back down, and a simple 3-tap kernel $[-1, 0, 1]$ — a discrete derivative that calculates slope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7edb3c4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def convolve1d(signal, kernel):\n",
    "    ksize = len(kernel)\n",
    "    pad = ksize // 2\n",
    "    padded = np.pad(signal, pad, mode='constant')\n",
    "    flipped = kernel[::-1]\n",
    "    out = np.zeros_like(signal, dtype=np.float64)\n",
    "    for n in range(len(signal)):\n",
    "        window = padded[n:n + ksize]\n",
    "        out[n] = np.dot(window, flipped)\n",
    "    return out\n",
    "\n",
    "signal = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1], dtype=np.float64)\n",
    "kernel = 0.5 * np.array([1, 0, -1], dtype=np.float64) # pre-flipped; scaling factor ensures that result is slope\n",
    "\n",
    "mine = convolve1d(signal, kernel)\n",
    "reference = np.convolve(signal, kernel, mode='same')\n",
    "print('signal:   ', signal)\n",
    "print('mine:     ', mine)\n",
    "print('np.convolve:', reference)\n",
    "print('match:', np.allclose(mine, reference))\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n",
    "axes[0].plot(signal, marker='o')\n",
    "axes[0].set_title('Signal')\n",
    "axes[1].plot(mine, marker='o', color='tab:orange')\n",
    "axes[1].axhline(0, color='gray', linewidth=0.8)\n",
    "axes[1].set_title('Signal * [1, 0, -1]')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50680868",
   "metadata": {},
   "source": [
    "Note the output is positive while the signal is rising, zero at the peak, and negative while it's falling — exactly what a derivative-like kernel should do."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86b6a4dd",
   "metadata": {},
   "source": [
    "## 2D convolution\n",
    "\n",
    "The 2D version of convolution (what actually gets used on images) is the same flip-and-slide idea as the 1D case above, just over two axes:\n",
    "\n",
    "$$(I * K)(x, y) = \\sum_{i}\\sum_{j} K(i, j)\\, I(x - i,\\, y - j)$$\n",
    "\n",
    "with image $I$ and kernel $K$.  Just as with the 1D version, the key detail is the *minus* signs: the kernel is flipped (rotated 180°) before it's slid across the image. This matters whenever the kernel is not symmetric.\n",
    "\n",
    "**Correlation** is the same idea without the flip:\n",
    "\n",
    "$$(I \\star K)(x, y) = \\sum_{i}\\sum_{j} K(i, j)\\, I(x + i,\\, y + j)$$\n",
    "\n",
    "In practice, most image-processing libraries — including OpenCV's `cv2.filter2D` — actually implement *correlation*, not convolution, even though people casually call it \"convolving with a kernel.\" For symmetric kernels (box blur, Gaussian) the two are identical, so the distinction rarely matters in practice. But it's worth knowing the difference exists."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88eb2dc2",
   "metadata": {},
   "source": [
    "## Convolution from scratch\n",
    "\n",
    "To implement 2D convolution: flip the kernel both horizontally and vertically, pad the image so the kernel never runs off the edge, then at every position multiply the flipped kernel elementwise by the pixels underneath it and sum. The code below implements this directly — except that for speed it loops over the kernel's (few) entries rather than the image's (many) pixels, shifting and accumulating the whole padded image at once for each kernel weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efdfecb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def convolve2d(image, kernel, border=cv2.BORDER_REFLECT101):\n",
    "    \"\"\"True convolution (kernel is flipped), single-channel, float output.\"\"\"\n",
    "    kh, kw = kernel.shape\n",
    "    pad_h, pad_w = kh // 2, kw // 2\n",
    "    flipped = kernel[::-1, ::-1]\n",
    "    padded = cv2.copyMakeBorder(image, pad_h, pad_h, pad_w, pad_w, border)\n",
    "\n",
    "    out = np.zeros(image.shape, dtype=np.float64)\n",
    "    for i in range(kh):  # loop over kernel values for speed (kernel is smaller than image)\n",
    "        for j in range(kw):\n",
    "            out += flipped[i, j] * padded[i:i + image.shape[0], j:j + image.shape[1]].astype(np.float64)\n",
    "    return out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57abeb5c",
   "metadata": {},
   "source": [
    "### Sanity check against OpenCV\n",
    "\n",
    "Since `cv2.filter2D` computes correlation, we feed it the *pre-flipped* kernel to get a true convolution result to compare against."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dbffdfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_img = np.zeros((20, 20), dtype=np.uint8)\n",
    "test_img[3:17, 3:6] = 200   # a small \"corner\" shape (no symmetry), so flipping is visible\n",
    "test_img[3:6, 3:14] = 200\n",
    "\n",
    "asymmetric_kernel = np.array([[1, 2, -1],\n",
    "                               [0, 1,  3],\n",
    "                               [-2, 1, 0]], dtype=np.float64)\n",
    "\n",
    "mine = convolve2d(test_img, asymmetric_kernel)\n",
    "reference = cv2.filter2D(test_img, cv2.CV_64F, asymmetric_kernel[::-1, ::-1],\n",
    "                          borderType=cv2.BORDER_REFLECT101)\n",
    "\n",
    "print('max abs difference:', np.abs(mine - reference).max())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c795c3d1",
   "metadata": {},
   "source": [
    "### Convolution vs. correlation, made visible\n",
    "\n",
    "For a symmetric kernel, convolving and correlating give *identical* results. For an asymmetric kernel, the results are different — but most asymmetric kernels in practice cause a simple sign flip."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebbe2d62",
   "metadata": {},
   "outputs": [],
   "source": [
    "correlated = cv2.filter2D(test_img, cv2.CV_64F, asymmetric_kernel, borderType=cv2.BORDER_REFLECT101)\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(9, 3.5))\n",
    "for ax, im, title in zip(axes, [test_img, mine, correlated],\n",
    "                          ['Original', 'Convolution\\n(kernel flipped)', 'Correlation\\n(kernel not flipped)']):\n",
    "    ax.imshow(im, cmap='gray')\n",
    "    ax.set_title(title, fontsize=9)\n",
    "    ax.axis('off')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print('convolution == correlation with kernel rotated 180?',\n",
    "      np.allclose(mine, cv2.filter2D(test_img, cv2.CV_64F, asymmetric_kernel[::-1, ::-1],\n",
    "                                      borderType=cv2.BORDER_REFLECT101)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05b4428e",
   "metadata": {},
   "source": [
    "## Boundary handling\n",
    "\n",
    "Near the edges, the kernel hangs off the image. `convolve2d` uses `cv2.copyMakeBorder` to pad first; the padding mode changes the result at the border. Common choices: constant (zero) padding, replicate (extend the edge pixel), and reflect (mirror across the edge, `REFLECT101` avoids duplicating the edge pixel itself). In this example, the inner $4 \\times 4$ array of each output is identical to the original array — only the outer ring is different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e78f3fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "small = np.arange(1, 17, dtype=np.float64).reshape(4, 4)\n",
    "\n",
    "for name, mode in [('CONSTANT (zero)', cv2.BORDER_CONSTANT),\n",
    "                    ('REPLICATE', cv2.BORDER_REPLICATE),\n",
    "                    ('REFLECT101', cv2.BORDER_REFLECT101)]:\n",
    "    padded = cv2.copyMakeBorder(small, 1, 1, 1, 1, mode)\n",
    "    print(f'{name}:\\n{padded}\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06ee81e6",
   "metadata": {},
   "source": [
    "## Classic filters as kernels\n",
    "\n",
    "Once we can convolve, most filters are just \"pick a kernel\":\n",
    "\n",
    "- **Box blur**: every neighbor weighted equally — averages out noise but blurs edges.\n",
    "- **Gaussian blur**: neighbors weighted by a bell curve — smoother falloff, less \"boxy\" artifacting than a box blur.\n",
    "- **Sharpen**: boosts the center pixel relative to its neighbors.\n",
    "- **Sobel (edge detection)**: an asymmetric kernel that responds strongly to intensity changes in one direction — this is a case where the convolution-vs-correlation flip actually matters, since the kernel is asymmetric (results in a sign flip)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f92b873e",
   "metadata": {},
   "outputs": [],
   "source": [
    "photo_like = np.zeros((100, 100), dtype=np.uint8)\n",
    "cv2.rectangle(photo_like, (20, 20), (80, 80), 200, -1)\n",
    "cv2.circle(photo_like, (50, 50), 20, 60, -1)\n",
    "rng = np.random.default_rng(0)\n",
    "photo_like = photo_like.astype(np.float64) + rng.normal(0, 12, photo_like.shape)\n",
    "photo_like = np.clip(photo_like, 0, 255)\n",
    "\n",
    "box = np.ones((5, 5)) / 25\n",
    "\n",
    "gx, gy = np.meshgrid(np.arange(5) - 2, np.arange(5) - 2)\n",
    "sigma = 1.0\n",
    "gaussian = np.exp(-(gx**2 + gy**2) / (2 * sigma**2))\n",
    "gaussian /= gaussian.sum()\n",
    "\n",
    "sharpen = np.array([[0, -1, 0],\n",
    "                     [-1, 5, -1],\n",
    "                     [0, -1, 0]], dtype=np.float64)\n",
    "\n",
    "sobel_x = np.array([[-1, 0, 1],\n",
    "                     [-2, 0, 2],\n",
    "                     [-1, 0, 1]], dtype=np.float64)\n",
    "\n",
    "results = {\n",
    "    'Original (noisy)': photo_like,\n",
    "    'Box blur': convolve2d(photo_like, box),\n",
    "    'Gaussian blur': convolve2d(photo_like, gaussian),\n",
    "    'Sharpen': convolve2d(photo_like, sharpen),\n",
    "    'Sobel x (edges)': convolve2d(photo_like, sobel_x),\n",
    "}\n",
    "\n",
    "fig, axes = plt.subplots(1, 5, figsize=(15, 3.5))\n",
    "for ax, (title, im) in zip(axes, results.items()):\n",
    "    ax.imshow(im, cmap='gray')\n",
    "    ax.set_title(title, fontsize=9)\n",
    "    ax.axis('off')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34da84e6",
   "metadata": {},
   "source": [
    "## Separable kernels: a speed trick\n",
    "\n",
    "The Gaussian kernel above is **separable**: it can be written as the outer product of two 1D kernels, $K = k_x \\, k_y^\\top$. That means convolving with the full $n \\times n$ kernel is equivalent to convolving with a $1 \\times n$ kernel, then an $n \\times 1$ kernel — turning $O(n^2)$ work per pixel into $O(2n)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f66c9875",
   "metadata": {},
   "outputs": [],
   "source": [
    "k1d = np.exp(-(np.arange(5) - 2) ** 2 / (2 * sigma**2))\n",
    "k1d /= k1d.sum()\n",
    "\n",
    "outer_product = np.outer(k1d, k1d)\n",
    "print('2D kernel == outer product of 1D kernels?', np.allclose(gaussian, outer_product))\n",
    "\n",
    "separable_result = convolve2d(convolve2d(photo_like, k1d.reshape(1, -1)), k1d.reshape(-1, 1))\n",
    "full_result = convolve2d(photo_like, gaussian)\n",
    "\n",
    "print('separable == full 2D convolution?', np.allclose(separable_result, full_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dc67665",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "1. Verify that the box kernel is also separable, by writing it as an outer product of two 1D uniform kernels.\n",
    "2. `sobel_x` above is not symmetric. Compute both the convolution and correlation of `photo_like` with it, and describe how the two outputs differ (hint: think about which direction of intensity change each responds to).\n",
    "3. Time `convolve2d` against `cv2.filter2D` on a 300x300 image with a 15x15 Gaussian kernel using `%timeit`. Then time the separable two-pass version against the full 2D version. How much faster is each speedup?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.x"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}