{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c8b8400b",
   "metadata": {},
   "source": [
    "# Lesson 7: Distance Measures\n",
    "\n",
    "\"Distance\" on a pixel grid is not as simple as it sounds. We look at three different flavors:\n",
    "\n",
    "1. **Point-to-point distance metrics** (Manhattan, chessboard, and Euclidean) — which measure how far apart two pixels are.\n",
    "2. **Curve-length estimators** (Freeman and a refined weighted estimator) — which measure the length of a *digitized curve* (e.g. the boundary of a blob) from its chain code.\n",
    "3. **Chamfer distance transform** — a fast, image-wide approximation to the Euclidean distance transform, computed with two raster-scan passes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48f25c31",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ad75fb2",
   "metadata": {},
   "source": [
    "## Part 1: Point-to-point distance metrics\n",
    "\n",
    "Given two pixels $p=(x_1,y_1)$ and $q=(x_2,y_2)$, let $d_x = x_2-x_1$ and $d_y=y_2-y_1$. Three common metrics:\n",
    "\n",
    "- **Manhattan / city-block ($L_1$)**: distance if you can move only along the grid axes (like walking city blocks).\n",
    "\n",
    "    $D_4(p,q) = |d_x| + |d_y|$ \n",
    "- **Chessboard ($L_\\infty$)**: number of king moves on a chessboard (diagonal steps are \"free\").\n",
    "\n",
    "    $D_8(p,q) = \\max(|d_x|,\\,|d_y|)$ \n",
    "- **Euclidean ($L_2$)**: ordinary straight-line distance.\n",
    "\n",
    "    $D_E(p,q) = \\sqrt{d_x^2+d_y^2}$\n",
    "\n",
    "The names $D_4$ and $D_8$ come from the fact that each is the shortest path length when only 4-connected (axis) or 8-connected (axis + diagonal) moves are allowed, one unit per move."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3baf8215",
   "metadata": {},
   "outputs": [],
   "source": [
    "def manhattan(dx, dy):\n",
    "    return np.abs(dx) + np.abs(dy)\n",
    "\n",
    "def chessboard(dx, dy):\n",
    "    return np.maximum(np.abs(dx), np.abs(dy))\n",
    "\n",
    "def euclidean(dx, dy):\n",
    "    return np.sqrt(dx**2 + dy**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ca3505d",
   "metadata": {},
   "source": [
    "### Visualizing the three metrics as distance fields\n",
    "\n",
    "For every pixel in a grid, we compute its distance to the center pixel under each metric and display it as an image. The metrics agree only along the axes — everywhere else they diverge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "923c6bb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 101\n",
    "half = size // 2\n",
    "yy, xx = np.mgrid[-half:half + 1, -half:half + 1]\n",
    "\n",
    "fields = {\n",
    "    'Manhattan ($D_4$)': manhattan(xx, yy),\n",
    "    'Chessboard ($D_8$)': chessboard(xx, yy),\n",
    "    'Euclidean': euclidean(xx, yy),\n",
    "}\n",
    "\n",
    "# Compute a single vmin/vmax to prevent each imshow() from auto-scaling \n",
    "# (Manhattan's max is 100, Euclidean's is ~71, chessboard's is 50).\n",
    "vmax = max(field.max() for field in fields.values())\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 4))\n",
    "for ax, (name, field) in zip(axes, fields.items()):\n",
    "    im = ax.imshow(field, cmap='viridis', vmin=0, vmax=vmax)\n",
    "    ax.contour(field, levels=8, colors='white', linewidths=0.5)\n",
    "    ax.set_title(name)\n",
    "    ax.axis('off')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28f22adb",
   "metadata": {},
   "source": [
    "The contour lines reveal each metric's characteristic shape: Manhattan distance forms diamonds, chessboard distance forms squares, and Euclidean distance forms circles — the one shape most people intuitively think of as \"equal distance away\"."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa8c0fa3",
   "metadata": {},
   "source": [
    "## Part 2: Estimating curve length from a chain code\n",
    "\n",
    "Now suppose we want to estimate the length of a *curve*, not just the distance between two points. Let's create a test image of a circle with known radius ($60$) and known circumference ($120 \\pi = 377$). Note that the number of boundary pixels ($336$) is a poor estimate of the circumference, since $336 \\neq 377$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65de2555",
   "metadata": {},
   "outputs": [],
   "source": [
    "radius = 60\n",
    "img = np.zeros((150, 150), dtype=np.uint8)\n",
    "cv2.circle(img, (75, 75), radius, 255, 1)\n",
    "\n",
    "contours, _ = cv2.findContours(img, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "pts = contours[0][:, 0, :]\n",
    "\n",
    "plt.imshow(img, cmap='gray')\n",
    "plt.title(f'Digitized circle, radius={radius}, {len(pts)} boundary pixels')\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa484384",
   "metadata": {},
   "source": [
    "A **Freeman chain code** describes a digital curve as a sequence of unit steps in 8 possible directions (0–7, 45° apart). Four of the eight directions are *even* (horizontal/vertical, true length 1) and four are *odd* (diagonal, true length $\\sqrt{2}$).\n",
    "\n",
    "If a curve's chain code has $N_e$ even steps and $N_o$ odd steps, **Freeman's estimate** of its length is simply\n",
    "\n",
    "$$L_{\\text{Freeman}} = N_e + \\sqrt{2}\\,N_o$$\n",
    "\n",
    "This sounds right, but it systematically *overestimates* smooth curves: a real circle's boundary, when digitized, alternates through many short zig-zag \"staircase\" steps that add up to more than the true arc length.\n",
    "\n",
    "To fix this problem, the exact weights $(1, \\sqrt{2})$ can be replaced with weights $(a, b)$ tuned to minimize the average error over many digitized curves — an idea sometimes associated with Kimura's chain-code length correction. Commonly cited values (they vary slightly by source) are approximately $a \\approx 0.948$ and $b \\approx 1.343$:\n",
    "\n",
    "$$L_{\\text{corrected}} = a\\,N_e + b\\,N_o$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0335c88b",
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = np.diff(np.vstack([pts, pts[:1]]), axis=0)  # step from each boundary pixel to the next\n",
    "is_diagonal = (steps[:, 0] != 0) & (steps[:, 1] != 0)\n",
    "\n",
    "N_o = int(np.sum(is_diagonal))    # odd (diagonal) steps\n",
    "N_e = int(np.sum(~is_diagonal))   # even (axis-aligned) steps\n",
    "\n",
    "true_length = 2 * np.pi * radius\n",
    "L_freeman = N_e + np.sqrt(2) * N_o\n",
    "\n",
    "a, b = 0.948, 1.343\n",
    "L_corrected = a * N_e + b * N_o\n",
    "\n",
    "print(f'N_e (even/axis steps)     = {N_e}')\n",
    "print(f'N_o (odd/diagonal steps)  = {N_o}')\n",
    "print(f'number of boundary pixels = {len(pts):.1f}')\n",
    "print(f'true circumference        = {true_length:.1f}')\n",
    "print(f'Freeman estimate          = {L_freeman:.1f}  (error {100*(L_freeman-true_length)/true_length:+.1f}%)')\n",
    "print(f'Corrected estimate        = {L_corrected:.1f}  (error {100*(L_corrected-true_length)/true_length:+.1f}%)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5ff4aeb",
   "metadata": {},
   "source": [
    "The corrected weights bring the estimate much closer to the true circumference, without needing anything more than a count of even and odd steps — useful when you want a cheap length estimate directly from a chain code."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b79e078f",
   "metadata": {},
   "source": [
    "## Part 3: Chamfer distance transform\n",
    "\n",
    "A **distance transform** labels every pixel of a binary image with its distance to the nearest foreground (or background) pixel — useful for skeletonization, shape matching, and path planning. Computing it *exactly* with Euclidean distance requires comparing every pixel against every foreground pixel (slow), or a more careful algorithm.\n",
    "\n",
    "The **chamfer distance transform** is a fast approximation: instead of a single global search, it propagates local distances across the image in just two raster-scan passes (top-left to bottom-right, then bottom-right to top-left), adding a small fixed cost for each step to a neighbor. Using integer weights $a=3$ for axis-aligned neighbors and $b=4$ for diagonal neighbors (then dividing by 3 at the end) gives the classic **3-4 chamfer distance**, a good, cheap approximation to Euclidean distance — the same idea as the Freeman/corrected chain-code weights from Part 2, but applied to every pixel instead of just a boundary.\n",
    "\n",
    "OpenCV's `cv2.distanceTransform` with `DIST_L2` computes the chamfer distance. The `maskSize` governs the quality of the approximation: a larger mask is more accurate but requires more compute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4a36798",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = np.zeros((100, 100), dtype=np.uint8)\n",
    "cv2.circle(img, (50, 50), 10, 255, -1)\n",
    "img = 255 - img  # we will compute distance from each background pixel to the nearest foreground (nonzero) pixel\n",
    "\n",
    "im_dist3 = cv2.distanceTransform(img, cv2.DIST_L2, maskSize = 3) # (3,4) chamfer algorithm \n",
    "im_dist5 = cv2.distanceTransform(img, cv2.DIST_L2, maskSize = 5) # more accurate algorithm\n",
    "\n",
    "error = im_dist3 - im_dist5\n",
    "print(f'max error  = {np.abs(error).max():.3f} pixels')\n",
    "print(f'mean error = {np.abs(error).mean():.3f} pixels')\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 4))\n",
    "for ax, field, title in zip(\n",
    "    axes,\n",
    "    [im_dist3, im_dist5, error],\n",
    "    ['Chamfer (3-4) approx.', 'More accurate estimate', 'Difference'],\n",
    "):\n",
    "    im = ax.imshow(field, cmap='viridis' if title != 'Difference' else 'coolwarm')\n",
    "    ax.set_title(title, fontsize=10)\n",
    "    ax.axis('off')\n",
    "    plt.colorbar(im, ax=ax, fraction=0.046)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "718dc8b6",
   "metadata": {},
   "source": [
    "The chamfer transform approximates the exact Euclidean distance transform, computed in two fast linear passes instead of an expensive nearest-neighbor search."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e63a8034",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "1. Repeat the circle experiment for a few different radii (e.g. 10, 20, 40, 80). Does the Freeman estimate's percent error stay roughly constant, or does it shrink as the circle gets bigger? Why might that be?\n",
    "2. The chessboard distance $D_8$ is the exact shortest path length when diagonal moves are allowed and cost 1. Can you construct a modified chessboard-like metric, using a diagonal cost of $\\sqrt{2}$ instead of 1, that behaves like a coarse local version of the chain-code length estimators above? Compare it to true Euclidean distance for a few points."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.x"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}