{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "28a970bc",
   "metadata": {},
   "source": [
    "# Lesson 6: Eigenvectors and Eigenvalues of the Moment Covariance Matrix\n",
    "\n",
    "In Lesson 5 we used central moments to find a blob's orientation and, with a bit of extra algebra, its semi-axis lengths. This lesson makes that connection explicit and general: the central moments of *any* blob (not just an ellipse) define a $2\\times2$ **covariance matrix**, and the eigenvectors/eigenvalues of that matrix directly give the orientation and size of the *equivalent ellipse* — the ellipse with the same area and same second-moment spread as the blob."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83313b1d",
   "metadata": {},
   "source": [
    "## Refresher: what are eigenvectors and eigenvalues?\n",
    "\n",
    "For a square matrix $A$, a nonzero vector $v$ is an **eigenvector** with **eigenvalue** $\\lambda$ if\n",
    "\n",
    "$$Av = \\lambda v$$\n",
    "\n",
    "In words: applying $A$ to $v$ doesn't rotate $v$ off its own line — it only *scales* it, by a factor of $\\lambda$. Most vectors get both rotated and scaled by $A$; eigenvectors are the special directions that only get scaled.\n",
    "\n",
    "For a $2\\times2$ **symmetric** matrix (like the covariance matrices we'll build below), something even nicer is guaranteed: there are always two eigenvectors, they are perpendicular to each other, and their eigenvalues are real numbers. Geometrically, $A$ takes a circle of unit vectors and stretches it into an ellipse whose axes point along the eigenvectors, with each semi-axis length equal to the corresponding eigenvalue. That picture — eigenvectors as axis *directions*, eigenvalues as axis *lengths* — is exactly what we'll exploit for shape analysis below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "027ae1ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "A = np.array([[3.0, 1.0],\n",
    "              [1.0, 1.5]])\n",
    "\n",
    "eigvals, eigvecs = np.linalg.eigh(A)  # ascending order\n",
    "print('eigenvalues :', eigvals)\n",
    "print('eigenvectors (columns):\\n', eigvecs)\n",
    "\n",
    "theta = np.linspace(0, 2 * np.pi, 200)\n",
    "circle = np.stack([np.cos(theta), np.sin(theta)])  # unit circle, as column vectors\n",
    "ellipse = A @ circle  # apply A to every point on the circle\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(4.5, 4.5))\n",
    "ax.plot(circle[0], circle[1], '--', color='gray', label='unit circle')\n",
    "ax.plot(ellipse[0], ellipse[1], color='#3498db', label='A @ circle')\n",
    "for val, vec in zip(eigvals, eigvecs.T):\n",
    "    ax.plot([0, val * vec[0]], [0, val * vec[1]], color='red', linewidth=2)\n",
    "ax.scatter(0, 0, color='black', zorder=5)\n",
    "ax.set_aspect('equal')\n",
    "ax.legend(loc='upper left', fontsize=8)\n",
    "ax.set_title('A circle stretched by A; red lines = eigenvectors x eigenvalues')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c07574e",
   "metadata": {},
   "source": [
    "### From one eigenvector to both at once: diagonalization\n",
    "\n",
    "$Av = \\lambda v$ above holds for a single eigenvector. Stack *both* eigenvectors as the columns of a matrix $P = \\begin{bmatrix}v_1 & v_2\\end{bmatrix}$, and both eigenvalues on the diagonal of $\\Lambda = \\begin{bmatrix}\\lambda_1 & 0\\\\0 & \\lambda_2\\end{bmatrix}$. Then $Av_1=\\lambda_1 v_1$ and $Av_2=\\lambda_2 v_2$, side by side, become one matrix equation:\n",
    "\n",
    "$$AP = P\\Lambda$$\n",
    "\n",
    "Because $A$ is symmetric, its eigenvectors are not just perpendicular but *orthonormal*, which makes $P$ an **orthogonal** matrix (inverting it is just a transpose: $P^{-1} = P^\\top$). That lets us solve for $A$ itself:\n",
    "\n",
    "$$A = P\\Lambda P^\\top$$\n",
    "\n",
    "This is the **diagonalization** of $A$, and geometrically it's three steps: $P^\\top$ *rotates* coordinates into the eigenvector frame, $\\Lambda$ *scales* along those (now axis-aligned) directions independently, and $P$ *rotates back*. Read in that rotated frame, $A$ is genuinely diagonal — the off-diagonal entry, which mixes $x$ and $y$ together, is exactly zero. When $A$ is a covariance matrix, that off-diagonal entry is $\\mu_{11}$: diagonalizing the covariance matrix and finding the ellipse's own natural axes are the same operation — both mean rotating into the one coordinate frame where the blob's spread in $x$ and $y$ no longer mixes at all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb1dd6c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "P = eigvecs\n",
    "Lam = np.diag(eigvals)\n",
    "\n",
    "print('A @ P:\\n', np.round(A @ P, 4))\n",
    "print('P @ Lambda:\\n', np.round(P @ Lam, 4))\n",
    "print('max |A@P - P@Lambda| (should be ~0):', np.abs(A @ P - P @ Lam).max())\n",
    "print()\n",
    "\n",
    "reconstructed = P @ Lam @ P.T\n",
    "print('P @ Lambda @ P.T (reconstructed A):\\n', np.round(reconstructed, 4))\n",
    "print('original A:\\n', A)\n",
    "print('max |reconstructed - A|:', np.abs(reconstructed - A).max())\n",
    "print()\n",
    "\n",
    "diagonalized = P.T @ A @ P  # A, read in the eigenvector frame\n",
    "print('P.T @ A @ P (A in its own eigenbasis -- should be diagonal):\\n', np.round(diagonalized, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06080955",
   "metadata": {},
   "source": [
    "## From central moments to a covariance matrix\n",
    "\n",
    "Treat a blob's pixels as samples from a 2D distribution. Its covariance matrix, in terms of the central moments $\\mu_{ij}$ and area $\\mu_{00}=m_{00}$, is\n",
    "\n",
    "$$A = \\frac{1}{m_{00}}\\begin{bmatrix}\\mu_{20} & \\mu_{11} \\\\ \\mu_{11} & \\mu_{02}\\end{bmatrix}$$\n",
    "\n",
    "This is exactly the same formula used for the covariance matrix of a scatter of points $(x,y)$, just weighted by pixel membership instead of by sample index. The off-diagonal $\\mu_{11}$ term is precisely the entry that mixing $x$ and $y$ leaves behind, so this is exactly the matrix we want to **diagonalize**: its eigenvectors will point along the blob's own major and minor axes — the directions of greatest and least spread — and its eigenvalues will measure how much spread there is along each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02fde223",
   "metadata": {},
   "outputs": [],
   "source": [
    "def covariance_from_moments(binary):\n",
    "    m = cv2.moments(binary, binaryImage=True)\n",
    "    cx, cy = m['m10'] / m['m00'], m['m01'] / m['m00']\n",
    "    cov = np.array([[m['mu20'], m['mu11']],\n",
    "                     [m['mu11'], m['mu02']]]) / m['m00']\n",
    "    return cov, (cx, cy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "affa28a2",
   "metadata": {},
   "source": [
    "## Eigen decomposition gives orientation and size\n",
    "\n",
    "`np.linalg.eigh` returns eigenvalues in ascending order along with their eigenvectors (as columns). For a filled ellipse with semi-axes $a \\ge b$, the eigenvalues work out to $\\lambda_{\\max} = a^2/4$ and $\\lambda_{\\min}=b^2/4$, so\n",
    "\n",
    "$$a = 2\\sqrt{\\lambda_{\\max}}, \\qquad b = 2\\sqrt{\\lambda_{\\min}}$$\n",
    "\n",
    "and the corresponding eigenvectors point along the major and minor axes.  Note that the angle of the major axis is the same as the angle computed in Lesson 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b688534",
   "metadata": {},
   "outputs": [],
   "source": [
    "def principal_axes(binary):\n",
    "    cov, center = covariance_from_moments(binary)\n",
    "    eigvals, eigvecs = np.linalg.eigh(cov)  # ascending order\n",
    "    semi_axes = 2 * np.sqrt(np.clip(eigvals, 0, None))\n",
    "    # reorder so index 0 is major (largest), index 1 is minor\n",
    "    order = [1, 0]\n",
    "    # eigvecs[:, order] has the two eigenvectors as its COLUMNS; transpose so that\n",
    "    # unpacking it below (which iterates over ROWS) hands back the two actual eigenvectors,\n",
    "    # not a scrambled mix of their x- and y-components\n",
    "    return center, semi_axes[order], eigvecs[:, order].T\n",
    "\n",
    "\n",
    "def draw_axes(ax, binary, color='red'):\n",
    "    center, (a, b), (v_major, v_minor) = principal_axes(binary)\n",
    "    cx, cy = center\n",
    "    ax.imshow(binary, cmap='gray')\n",
    "    for length, vec, lw in [(a, v_major, 2.5), (b, v_minor, 1.5)]:\n",
    "        dx, dy = length * vec\n",
    "        ax.plot([cx - dx, cx + dx], [cy - dy, cy + dy], c=color, linewidth=lw)\n",
    "    ax.scatter(cx, cy, c=color, marker='x', s=60)\n",
    "    ax.axis('off')\n",
    "    return center, (a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed793576",
   "metadata": {},
   "source": [
    "### Sanity check: a known ellipse\n",
    "\n",
    "As in Lesson 5, we draw an ellipse with known parameters and confirm the eigen-based estimate recovers them — but this time we plot *both* axes, not just the major one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "542d925b",
   "metadata": {},
   "outputs": [],
   "source": [
    "binary = np.zeros((200, 200), dtype=np.uint8)\n",
    "cv2.ellipse(binary, (100, 100), (70, 25), 30, 0, 360, 255, -1)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(4, 4))\n",
    "center, (a, b) = draw_axes(ax, binary)\n",
    "ax.set_title('Major (thick) and minor (thin) axes from eigenvectors')\n",
    "plt.show()\n",
    "\n",
    "print(f'recovered semi-axes: a={a:.1f}, b={b:.1f}  (drawn with 70, 25)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd541959",
   "metadata": {},
   "source": [
    "## It works on arbitrary shapes, too\n",
    "\n",
    "The real power of this approach is that it doesn't require the blob to be an ellipse at all. Every binary shape has *some* equivalent ellipse — the one that matches its area, centroid, and second-moment spread. This gives a compact 5-number summary (center, two semi-axis lengths, orientation) of an arbitrarily shaped blob."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89fba3ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "im_glasses = cv2.imread('../img/glasses_outline.png', cv2.IMREAD_GRAYSCALE)\n",
    "\n",
    "shapes = {'Glasses': im_glasses}\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 4))\n",
    "center, (a, b) = draw_axes(ax, im_glasses)\n",
    "ax.set_title(f'Glasses with axes\\na={a:.0f}, b={b:.0f}', fontsize=10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fdbb6b3",
   "metadata": {},
   "source": [
    "Notice the equivalent ellipse doesn't try to trace the shape's boundary — it summarizes the *distribution of mass* around the centroid, so it cuts across the narrow bridge between the two lenses rather than following it, and its outline falls outside the lenses near the top and bottom where the glasses taper."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2bf20e2",
   "metadata": {},
   "source": [
    "## Eccentricity: how elongated is a shape?\n",
    "\n",
    "The ratio of the eigenvalues (or semi-axes) gives a single number describing elongation, independent of orientation and overall size:\n",
    "\n",
    "$$\\text{eccentricity} = \\sqrt{1 - \\frac{\\lambda_{\\min}}{\\lambda_{\\max}}}$$\n",
    "\n",
    "This ranges from 0 (a circle, both axes equal) to nearly 1 (a very thin, elongated shape)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b063c9fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'{\"shape\":>10} {\"a\":>6} {\"b\":>6} {\"eccentricity\":>13}')\n",
    "for name, img in shapes.items():\n",
    "    _, (a, b), _ = principal_axes(img)\n",
    "    ecc = np.sqrt(1 - (b / a) ** 2)\n",
    "    print(f'{name:>10} {a:6.1f} {b:6.1f} {ecc:13.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "314a36d5",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "1. Draw a perfect circle and confirm its eccentricity is (close to) 0.\n",
    "2. Overlay the *actual* fitted ellipse outline (not just the axis lines) using `cv2.ellipse` with the center, `(2a, 2b)` as the full axes lengths, and the orientation angle from the eigenvectors. Compare it to `cv2.fitEllipse` applied to the shape's contour — do they agree?\n",
    "3. Run `principal_axes` on the filled blob of `glasses` (after flood fill)). How much do `a`, `b`, and the eccentricity change? Does that match your intuition for how moments depend on *where* mass sits, not just the overall silhouette?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.x"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}