{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9aa037fa",
   "metadata": {},
   "source": [
    "# Lesson 23: Projective Geometry\n",
    "\n",
    "Lesson 8 built a hierarchy of 2D transforms — Euclidean, similarity, affine — and noted that all three preserve parallel lines. This lesson covers the next, most general step: the **projective transform (homography)**, which models what a camera actually does when it looks at a flat surface from an angle, and which does *not* preserve parallelism. That's not a bug — it's exactly the phenomenon of a vanishing point, and it's the tool behind perspective correction and image stitching."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5716d2a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f140d2c",
   "metadata": {},
   "source": [
    "## Homogeneous coordinates\n",
    "\n",
    "Represent a 2D point $(x, y)$ as a 3-vector $(x, y, 1)$. Any $3\\times3$ matrix $H$ can then act on it by ordinary matrix multiplication; converting back to 2D means dividing by the third coordinate:\n",
    "\n",
    "$$\\begin{bmatrix}x'\\\\y'\\\\w'\\end{bmatrix} = H\\begin{bmatrix}x\\\\y\\\\1\\end{bmatrix}, \\qquad (x_{\\text{2D}}', y_{\\text{2D}}') = \\left(\\frac{x'}{w'}, \\frac{y'}{w'}\\right)$$\n",
    "\n",
    "Two big payoffs: first, a $3\\times3$ matrix can represent translation too (impossible with a $2\\times2$ matrix alone, which is why Lesson 8 needed a separate translation vector $t$). Second, that division by $w'$ is what lets a homography model genuine *perspective* effects — not just the affine transforms of Lesson 8, but photos where parallel lines converge toward a vanishing point, exactly what a real camera does when it looks at a flat surface from an angle. Also note that $H$ and $cH$ (any nonzero scalar multiple) represent the *exact same* transform, since the division cancels the scale — a homography has only 8 independent degrees of freedom, not 9."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "682b088c",
   "metadata": {},
   "source": [
    "## Parallel lines stop being parallel\n",
    "\n",
    "The updated hierarchy, extending Lesson 8's table:\n",
    "\n",
    "| Transform | Matrix form | Preserves |\n",
    "|---|---|---|\n",
    "| Affine | $\\begin{bmatrix}a&b&t_x\\\\c&d&t_y\\\\0&0&1\\end{bmatrix}$ | parallelism |\n",
    "| **Projective** | $\\begin{bmatrix}a&b&t_x\\\\c&d&t_y\\\\g&h&1\\end{bmatrix}$ | straight lines only |\n",
    "\n",
    "The only difference is a nonzero bottom row $(g, h)$ — and that alone is enough to destroy parallelism. We demonstrate directly: take two parallel horizontal lines and apply a homography with a nonzero bottom row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "000328d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "H = np.array([\n",
    "    [1,     0.2,   0],\n",
    "    [0.1,   1,     0],\n",
    "    [0.02,  0.015, 1],\n",
    "], dtype=np.float64)\n",
    "\n",
    "def apply_homography(points, H):\n",
    "    homogeneous = np.hstack([points, np.ones((len(points), 1))])\n",
    "    transformed = (H @ homogeneous.T).T\n",
    "    return transformed[:, :2] / transformed[:, 2:3]\n",
    "\n",
    "line1 = np.array([[0, 0], [1, 0]], dtype=np.float64)   # y = 0\n",
    "line2 = np.array([[0, 1], [1, 1]], dtype=np.float64)   # y = 1, parallel to line1\n",
    "\n",
    "l1_transformed = apply_homography(line1, H)\n",
    "l2_transformed = apply_homography(line2, H)\n",
    "\n",
    "print('line 1, transformed:', l1_transformed)\n",
    "print('line 2, transformed:', l2_transformed)\n",
    "\n",
    "def extended(p1, p2, length):\n",
    "    direction = (p2 - p1)\n",
    "    return p1 - length * direction, p2 + length * direction\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(9, 4))\n",
    "\n",
    "for line, color in [(line1, 'tab:blue'), (line2, 'tab:orange')]:\n",
    "    a, b = extended(line[0], line[1], length=3)\n",
    "    axes[0].plot([a[0], b[0]], [a[1], b[1]], color=color)\n",
    "    axes[0].plot(line[:, 0], line[:, 1], 'o', color=color)\n",
    "axes[0].set_xlim(-2, 4)\n",
    "axes[0].set_ylim(-2, 3)\n",
    "axes[0].set_title('Before: two parallel lines')\n",
    "\n",
    "for line, color in [(l1_transformed, 'tab:blue'), (l2_transformed, 'tab:orange')]:\n",
    "    a, b = extended(line[0], line[1], length=60)\n",
    "    axes[1].plot([a[0], b[0]], [a[1], b[1]], color=color)\n",
    "    axes[1].plot(line[:, 0], line[:, 1], 'o', color=color)\n",
    "axes[1].set_xlim(-5, 55)\n",
    "axes[1].set_ylim(-2, 8)\n",
    "axes[1].set_title('After: a homography with a nonzero bottom row')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce985328",
   "metadata": {},
   "source": [
    "## Vanishing points, computed two ways\n",
    "\n",
    "Where do these two lines actually meet? We can find it two ways: intersecting the transformed line segments directly, or — more elegantly — transforming the *point at infinity* in the original lines' shared direction $(1, 0)$, represented in homogeneous coordinates as $(1, 0, 0)$ (a nonzero third coordinate would make it a finite point; zero means \"infinitely far away\"). Applying $H$ to that point at infinity should land exactly on the vanishing point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5857483d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def line_intersection(p1, p2, p3, p4):\n",
    "    A = np.array([[p2[0] - p1[0], -(p4[0] - p3[0])],\n",
    "                  [p2[1] - p1[1], -(p4[1] - p3[1])]])\n",
    "    b = np.array([p3[0] - p1[0], p3[1] - p1[1]])\n",
    "    t = np.linalg.solve(A, b)[0]\n",
    "    return p1 + t * (p2 - p1)\n",
    "\n",
    "vanishing_point_direct = line_intersection(l1_transformed[0], l1_transformed[1],\n",
    "                                            l2_transformed[0], l2_transformed[1])\n",
    "\n",
    "point_at_infinity = np.array([1.0, 0.0, 0.0])  # direction (1,0), infinitely far away\n",
    "transformed_infinity = H @ point_at_infinity\n",
    "vanishing_point_via_infinity = transformed_infinity[:2] / transformed_infinity[2]\n",
    "\n",
    "print('vanishing point (line intersection): ', vanishing_point_direct)\n",
    "print('vanishing point (H @ infinity trick):', vanishing_point_via_infinity)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "376bed5d",
   "metadata": {},
   "source": [
    "## Lines in homogeneous coordinates\n",
    "\n",
    "Points aren't the only thing homogeneous coordinates represent elegantly — a line $ax+by+c=0$ is just as naturally the 3-vector $l=(a,b,c)$, with a point $p=(x,y,1)$ lying on it exactly when the inner product $l\\cdot p= l^\\top p = 0$. Two useful consequences follow from the cross product: the line through two points is $l = p_1\\times p_2$ (it's orthogonal to both, so both dot products vanish), and by the same logic in reverse, the intersection of two lines is $p = l_1\\times l_2$. `line_intersection` above solved a $2\\times2$ linear system by hand; the cross product does the same job in one line of code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a1bacac",
   "metadata": {},
   "outputs": [],
   "source": [
    "def line_through(p1, p2):\n",
    "    return np.cross([p1[0], p1[1], 1.0], [p2[0], p2[1], 1.0])\n",
    "\n",
    "def intersect(l1, l2):\n",
    "    p = np.cross(l1, l2)\n",
    "    return p[:2] / p[2]\n",
    "\n",
    "line1_h = line_through(l1_transformed[0], l1_transformed[1])\n",
    "line2_h = line_through(l2_transformed[0], l2_transformed[1])\n",
    "vanishing_point_cross = intersect(line1_h, line2_h)\n",
    "\n",
    "print('vanishing point (cross-product trick):', vanishing_point_cross)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edc0bafe",
   "metadata": {},
   "source": [
    "Matches the two earlier methods exactly. Points and lines (in the 2D plane) turn out to play completely symmetric roles in homogeneous coordinates — the same cross-product operation builds a line from two points or a point from two lines."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee34801f",
   "metadata": {},
   "source": [
    "## When does a homography actually apply?\n",
    "\n",
    "A homography exactly relates two photos in exactly two situations: photographing a flat surface from two different positions (the case used below), or photographing *any* 3D scene, flat or not, from the same camera position while only rotating the camera between shots. That second case is what makes Lesson 24's panorama stitching work on ordinary, non-planar scenes — a translating camera or a genuinely 3D scene viewed from two different positions has no single homography that relates the two images exactly."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8af1cd46",
   "metadata": {},
   "source": [
    "## Perspective correction: rectifying a photographed document\n",
    "\n",
    "The most common practical use of a homography: given 4 point correspondences (e.g. the 4 corners of a document, clicked by a user or found automatically), `cv2.getPerspectiveTransform` solves for the unique homography mapping one set of 4 points to the other, and `cv2.warpPerspective` applies it. The trapezoid shape a photographed rectangle takes on — edges that were parallel in real life visibly converging in the photo — is called **keystoning**, and it's exactly the vanishing-point effect above, now applied to a document instead of a pair of abstract lines. (Below, `H_rectify` is computed from the point correspondences in reverse order, but it's exactly `np.linalg.inv(H_distort)`, since undoing a homography and inverting its matrix are the same operation.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98d6d148",
   "metadata": {},
   "outputs": [],
   "source": [
    "document = np.full((300, 220, 3), 255, dtype=np.uint8)\n",
    "cv2.rectangle(document, (20, 20), (200, 280), (0, 0, 0), 3)\n",
    "cv2.putText(document, 'HELLO', (30, 150), cv2.FONT_HERSHEY_SIMPLEX, 1.2, (0, 0, 0), 3)\n",
    "\n",
    "h, w = document.shape[:2]\n",
    "corners = np.float32([[0, 0], [w, 0], [w, h], [0, h]])\n",
    "photographed_corners = np.float32([[45, 15], [175, 25], [210, 290], [5, 280]])  # strong keystone: top edge narrower than bottom\n",
    "\n",
    "H_distort = cv2.getPerspectiveTransform(corners, photographed_corners)\n",
    "photographed = cv2.warpPerspective(document, H_distort, (w, h))\n",
    "\n",
    "H_rectify = cv2.getPerspectiveTransform(photographed_corners, corners)\n",
    "rectified = cv2.warpPerspective(photographed, H_rectify, (w, h))\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(9, 4))\n",
    "for ax, im, title in zip(axes, [document, photographed, rectified],\n",
    "                          ['Original document', 'Photographed at an angle\\n(simulated)', 'Rectified\\n(from 4 corner clicks)']):\n",
    "    ax.imshow(cv2.cvtColor(im, cv2.COLOR_BGR2RGB))\n",
    "    ax.set_title(title, fontsize=9)\n",
    "    ax.axis('off')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "error = np.abs(rectified.astype(int) - document.astype(int))\n",
    "print(f'mean abs pixel difference after distort-then-rectify round trip: {error.mean():.2f}  (small resampling loss only)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22bc3037",
   "metadata": {},
   "source": "## Solving for a homography from many correspondences\n\n4 correspondences exactly determine a homography's 8 degrees of freedom, with no slack for error. With *more* than 4 (typically from automatic feature matching, Lesson 19), the problem becomes an overdetermined least-squares fit — the Direct Linear Transform (DLT) algorithm that `cv2.findHomography` implements, optionally wrapped in RANSAC (Lesson 22) to reject bad correspondences (mismatched features) as outliers. `cv2.findHomography` also normalizes the point coordinates internally before solving, for numerical conditioning — the same trick made explicit in Lesson 26's 8-point algorithm."
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f71a6b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(0)\n",
    "H_true = np.array([[1, 0.2, 10], [0.05, 1, 5], [0.0008, 0.0003, 1]])\n",
    "\n",
    "pts1 = rng.uniform(0, 200, (30, 2))\n",
    "pts2 = apply_homography(pts1, H_true)\n",
    "\n",
    "H_estimated, inlier_mask = cv2.findHomography(pts1, pts2, cv2.RANSAC, 3.0)\n",
    "\n",
    "print('true homography (scale-normalized):\\n', np.round(H_true / H_true[2, 2], 4))\n",
    "print('estimated homography:\\n', np.round(H_estimated, 4))\n",
    "print(f'inliers: {int(inlier_mask.sum())} / {len(inlier_mask)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d0e6b7",
   "metadata": {},
   "source": [
    "With clean correspondences, `findHomography` recovers `H_true` almost exactly — this is exactly the last step of a typical image-stitching pipeline: detect and match SIFT features between two overlapping photos (Lesson 19), then solve for the homography that aligns one onto the other."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "027c2137",
   "metadata": {},
   "source": [
    "## RANSAC in practice: homography estimation with bad matches\n",
    "\n",
    "Even though the code above passes `cv2.RANSAC` to `findHomography`, with clean correspondences there is nothing for it to reject — it makes no visible difference. What effect does RANSAC actually have? Since real feature matching (Lesson 19) inevitably produces mismatched features, even after a ratio test, we need a way to discard such outliers. To illustrate this, we inject 15 deliberately garbage correspondences alongside the 30 genuine ones and compare the RANSAC fit against a plain least-squares fit on the exact same contaminated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76115a5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "pts1_bad = rng.uniform(0, 200, (15, 2))\n",
    "pts2_bad = rng.uniform(0, 300, (15, 2))  # unrelated -- simulated bad matches\n",
    "\n",
    "pts1_contaminated = np.vstack([pts1, pts1_bad])\n",
    "pts2_contaminated = np.vstack([pts2, pts2_bad])\n",
    "\n",
    "H_ransac, mask = cv2.findHomography(pts1_contaminated, pts2_contaminated, cv2.RANSAC, 3.0)\n",
    "H_plain, _ = cv2.findHomography(pts1_contaminated, pts2_contaminated, 0)  # method=0: plain least squares, no outlier rejection\n",
    "\n",
    "def mean_reprojection_error(H, pts1, pts2):\n",
    "    return np.linalg.norm(apply_homography(pts1, H) - pts2, axis=1).mean()\n",
    "\n",
    "print(f'inliers found by RANSAC: {int(mask.sum())} / {len(mask)}  (30 correspondences were genuinely correct)')\n",
    "print()\n",
    "print(f'mean reprojection error on the TRUE inliers:')\n",
    "print(f'  RANSAC fit: {mean_reprojection_error(H_ransac, pts1, pts2):.2e} pixels')\n",
    "print(f'  plain fit:  {mean_reprojection_error(H_plain, pts1, pts2):.2e} pixels')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4e4ccb3",
   "metadata": {},
   "source": [
    "RANSAC finds exactly the 30 genuine correspondences and reconstructs the homography to essentially machine precision. The plain least-squares fit, given the exact same data, is off by many orders of magnitude more error — effectively useless for anything requiring pixel-level accuracy — because 15 bad matches out of 45 (33%) was more than enough to noticeably corrupt an unweighted sum-of-squares fit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d67ebd89",
   "metadata": {},
   "source": [
    "## Looking ahead: projective geometry in 3D\n",
    "\n",
    "Everything in this lesson has been 2D-to-2D: a homography relates two planes (images or flat surfaces in the world). Lesson 25 extends the same projective machinery one dimension further, deriving the camera's projection matrix that maps a 3D scene onto a 2D image in the first place."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ffc2ada",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "1. Add zero-mean Gaussian noise (e.g. std 2 pixels) to `pts2` before calling `cv2.findHomography`. How much does the estimated homography drift from `H_true`, and does increasing the number of correspondences (say, from 30 to 200) reduce that drift?\n",
    "2. Increase the number of injected bad matches (`pts1_bad`/`pts2_bad`) from 15 until RANSAC starts to fail. Using the iteration-count formula from Lesson 22 (minimal sample size 4 for a homography), roughly what outlier fraction is that, and how many iterations would it theoretically require to stay 99% confident of success?\n",
    "3. In the vanishing-point demo, change the two lines to be *vertical* instead of horizontal (e.g. `x=0` and `x=1`), and predict, then verify, where their vanishing point ends up using the point-at-infinity trick with direction `(0, 1, 0)`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.x"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}