Image source: BlendedMVS (CC BY 4.0)
" ] }, { "cell_type": "markdown", "id": "20a17371", "metadata": {}, "source": [ "## Step 1: feature matching\n", "\n", "SIFT + ratio test (Lesson 20), run directly on the two photos." ] }, { "cell_type": "code", "execution_count": null, "id": "323b120e", "metadata": {}, "outputs": [], "source": [ "gray1 = cv2.cvtColor(view1, cv2.COLOR_BGR2GRAY)\n", "gray2 = cv2.cvtColor(view2, cv2.COLOR_BGR2GRAY)\n", "\n", "sift = cv2.SIFT_create()\n", "kp1, des1 = sift.detectAndCompute(gray1, None)\n", "kp2, des2 = sift.detectAndCompute(gray2, None)\n", "\n", "bf = cv2.BFMatcher()\n", "raw_matches = bf.knnMatch(des1, des2, k=2)\n", "good_matches = [m for m, n in raw_matches if m.distance < 0.75 * n.distance]\n", "\n", "x1 = np.float32([kp1[m.queryIdx].pt for m in good_matches])\n", "x2 = np.float32([kp2[m.trainIdx].pt for m in good_matches])\n", "\n", "print(f'keypoints: {len(kp1)} (view 1), {len(kp2)} (view 2)')\n", "print(f'good matches after ratio test: {len(good_matches)}')\n", "\n", "match_vis = cv2.drawMatches(view1, kp1, view2, kp2, good_matches, None,\n", " flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)\n", "plt.figure(figsize=(10, 4))\n", "plt.imshow(cv2.cvtColor(match_vis, cv2.COLOR_BGR2RGB))\n", "plt.title('Feature matches')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "313721d7", "metadata": {}, "source": [ "## Step 2: robust essential matrix and pose recovery\n", "\n", "Feed the matched points through the same robust pipeline built in Lessons 23 and 27 to recover the essential matrix and the cameras' relative pose." ] }, { "cell_type": "code", "execution_count": null, "id": "f86ff230", "metadata": {}, "outputs": [], "source": [ "E, inlier_mask = cv2.findEssentialMat(x1, x2, K, method=cv2.RANSAC, threshold=1.0)\n", "_, R_estimated, t_estimated, _ = cv2.recoverPose(E, x1, x2, K)\n", "\n", "print(f'inliers: {int(inlier_mask.sum())} / {len(inlier_mask)}')\n", "print('recovered rotation:\\n', np.round(R_estimated, 4))\n", "print('true rotation:\\n', np.round(R_true, 4))\n", "print()\n", "print('recovered translation direction:', np.round(t_estimated.ravel(), 4))\n", "print('true translation direction: ', np.round(t_true / np.linalg.norm(t_true), 4))" ] }, { "cell_type": "markdown", "id": "3b8cb849", "metadata": {}, "source": [ "As in Lesson 28, `recoverPose` returns a **unit-length** translation direction — the actual baseline distance between the cameras is fundamentally unrecoverable from image correspondences alone. We'll come back to that." ] }, { "cell_type": "markdown", "id": "9f129208", "metadata": {}, "source": [ "## Step 3: triangulation\n", "\n", "Given two camera projection matrices $P_1, P_2$ (each $3\\times4$, mapping a 3D point to a 2D image point in homogeneous coordinates) and a matched pixel pair $(x_1, x_2)$, we want the 3D point $X$ satisfying both $x_1 \\propto P_1 X$ and $x_2 \\propto P_2 X$. Each view contributes 2 independent linear equations in $X$'s 4 homogeneous unknowns (cross-multiplying out the unknown scale factor), by the same DLT/SVD recipe used for homographies (Lesson 24) and the fundamental matrix (Lesson 28):\n", "\n", "$$A = \\begin{bmatrix} u_1 P_1^{(3)} - P_1^{(1)} \\\\ v_1 P_1^{(3)} - P_1^{(2)} \\\\ u_2 P_2^{(3)} - P_2^{(1)} \\\\ v_2 P_2^{(3)} - P_2^{(2)} \\end{bmatrix}, \\qquad AX = 0$$\n", "\n", "where $P^{(i)}$ denotes row $i$ of $P$. The null space of $A$ (the right singular vector associated with the smallest singular value) gives $X$." ] }, { "cell_type": "code", "execution_count": null, "id": "8aba34b9", "metadata": {}, "outputs": [], "source": [ "def triangulate_dlt(P1, P2, x1, x2):\n", " points = []\n", " for (u1, v1), (u2, v2) in zip(x1, x2):\n", " A = np.array([\n", " u1 * P1[2] - P1[0],\n", " v1 * P1[2] - P1[1],\n", " u2 * P2[2] - P2[0],\n", " v2 * P2[2] - P2[1],\n", " ])\n", " _, _, Vt = np.linalg.svd(A)\n", " X = Vt[-1]\n", " points.append(X[:3] / X[3])\n", " return np.array(points)\n", "\n", "# validate against the dataset's TRUE camera poses: triangulate the RANSAC inliers with them,\n", "# and check that reprojecting the resulting 3D points lands back on the original pixels\n", "inliers = inlier_mask.ravel().astype(bool)\n", "x1_in, x2_in = x1[inliers], x2[inliers]\n", "\n", "recon_true_cameras = triangulate_dlt(P1, P2_true, x1_in, x2_in)\n", "\n", "cv_points_4d = cv2.triangulatePoints(P1, P2_true, x1_in.T, x2_in.T)\n", "cv_points_3d = (cv_points_4d[:3] / cv_points_4d[3]).T\n", "print(f'max diff vs. cv2.triangulatePoints: {np.abs(recon_true_cameras - cv_points_3d).max():.2e}')\n", "\n", "# cheirality check: a correctly-matched 3D point must be in front of BOTH cameras --\n", "# a real pipeline always applies this, since epipolar-only RANSAC can still admit a bad match\n", "depth1 = recon_true_cameras[:, 2]\n", "depth2 = (R_true @ recon_true_cameras.T + t_true.reshape(3, 1))[2]\n", "valid = (depth1 > 0) & (depth2 > 0)\n", "x1_in, x2_in, recon_true_cameras = x1_in[valid], x2_in[valid], recon_true_cameras[valid]\n", "print(f'kept {valid.sum()} / {len(valid)} inliers after cheirality check')\n", "\n", "def reprojection_error(P, X, x):\n", " X_h = np.hstack([X, np.ones((len(X), 1))])\n", " proj = (P @ X_h.T).T\n", " proj = proj[:, :2] / proj[:, 2:3]\n", " return np.linalg.norm(proj - x, axis=1)\n", "\n", "err1 = reprojection_error(P1, recon_true_cameras, x1_in)\n", "err2 = reprojection_error(P2_true, recon_true_cameras, x2_in)\n", "all_err = np.concatenate([err1, err2])\n", "print(f'reprojection error using TRUE camera poses: mean {all_err.mean():.3f} px, max {all_err.max():.3f} px')" ] }, { "cell_type": "markdown", "id": "3882ef0f", "source": "### Cheirality: a free sanity check\n\nDecomposing $E$ into $(R, t)$ back in Step 2 actually yields four mathematically valid sign combinations; only one places the observed points in front of *both* cameras, and `cv2.recoverPose` uses exactly that fact — the **cheirality constraint** — to pick the right one internally, by triangulating a few points and checking their depth. The check above applies the same idea explicitly: the epipolar constraint alone doesn't rule out a correspondence that triangulates *behind* a camera, so a mismatched point can still pass RANSAC's inlier test. Requiring positive depth in both views is a free extra filter on top of that — it's what dropped the one inlier flagged above.", "metadata": {} }, { "cell_type": "markdown", "id": "3d4fc872", "metadata": {}, "source": [ "## Putting it together: reconstruction using *estimated* pose\n", "\n", "Now the real test: triangulate using the camera matrix built from `recoverPose`'s *estimated* $(R, t)$, not the dataset's true pose. Because $t$ was only recovered up to scale, so is the reconstruction — every 3D point comes out a fixed factor smaller than reality. Multiplying by the true baseline length (a bit of cheating, just for evaluation) should recover the correct metric scene." ] }, { "cell_type": "code", "execution_count": null, "id": "6db8e4fd", "metadata": {}, "outputs": [], "source": [ "P2_estimated = K @ np.hstack([R_estimated, t_estimated.reshape(3, 1)])\n", "reconstruction_unit_scale = triangulate_dlt(P1, P2_estimated, x1_in, x2_in)\n", "\n", "true_baseline = np.linalg.norm(t_true)\n", "reconstruction_metric = reconstruction_unit_scale * true_baseline\n", "\n", "error = np.linalg.norm(reconstruction_metric - recon_true_cameras, axis=1)\n", "print(f'true baseline: {true_baseline:.3f}')\n", "print(f'mean 3D reconstruction error (vs. triangulation from true poses): {error.mean():.3f}')\n", "print(f'max 3D reconstruction error: {error.max():.3f}')" ] }, { "cell_type": "markdown", "id": "fd6a74a3", "metadata": {}, "source": [ "With real (imperfect) correspondences, the entire pipeline — essential matrix, pose, triangulation — reconstructs the scene to within a small fraction of the camera baseline, using nothing but 2D pixel correspondences and known intrinsics, plus *one* external number (the true baseline) to fix the scale ambiguity. In practice, that scale reference might come from a known object size in the scene, a second sensor (GPS, IMU, LiDAR), or a calibrated stereo rig (Lesson 22) instead of two arbitrary independent cameras." ] }, { "cell_type": "markdown", "id": "df0ab549", "metadata": {}, "source": [ "### Visualizing the reconstruction\n", "\n", "For a denser, more recognizable point cloud than the strict RANSAC-inlier set used for the numeric check above, we loosen the feature-matching thresholds (`SIFT_create(contrastThreshold=0.01)`, ratio 0.85) to pull in more — slightly noisier — correspondences, then color each 3D point with its actual pixel color from view 1." ] }, { "cell_type": "code", "execution_count": null, "id": "42b6c294", "metadata": {}, "outputs": [], "source": [ "sift_dense = cv2.SIFT_create(contrastThreshold=0.01)\n", "kp1d, des1d = sift_dense.detectAndCompute(gray1, None)\n", "kp2d, des2d = sift_dense.detectAndCompute(gray2, None)\n", "raw_dense = bf.knnMatch(des1d, des2d, k=2)\n", "good_dense = [m for m, n in raw_dense if m.distance < 0.85 * n.distance]\n", "\n", "x1d = np.float32([kp1d[m.queryIdx].pt for m in good_dense])\n", "x2d = np.float32([kp2d[m.trainIdx].pt for m in good_dense])\n", "\n", "_, dense_mask = cv2.findEssentialMat(x1d, x2d, K, method=cv2.RANSAC, threshold=1.0)\n", "dense_inliers = dense_mask.ravel().astype(bool)\n", "x1d, x2d = x1d[dense_inliers], x2d[dense_inliers]\n", "\n", "recon_dense_true = triangulate_dlt(P1, P2_true, x1d, x2d)\n", "recon_dense_estimated = triangulate_dlt(P1, P2_estimated, x1d, x2d) * true_baseline\n", "\n", "# cheirality check again, same as for the sparse set above\n", "depth1 = recon_dense_true[:, 2]\n", "depth2 = (R_true @ recon_dense_true.T + t_true.reshape(3, 1))[2]\n", "valid = (depth1 > 0) & (depth2 > 0)\n", "x1d, recon_dense_true, recon_dense_estimated = x1d[valid], recon_dense_true[valid], recon_dense_estimated[valid]\n", "\n", "rgb1 = cv2.cvtColor(view1, cv2.COLOR_BGR2RGB)\n", "px = np.clip(x1d[:, 0].astype(int), 0, rgb1.shape[1] - 1)\n", "py = np.clip(x1d[:, 1].astype(int), 0, rgb1.shape[0] - 1)\n", "point_colors = rgb1[py, px] / 255.0\n", "\n", "print(f'{len(x1d)} points in the denser cloud (vs. {len(x1_in)} used for the numeric check above)')" ] }, { "cell_type": "code", "execution_count": null, "id": "e9360d20", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(7, 6))\n", "ax = fig.add_subplot(111, projection='3d')\n", "ax.scatter(*recon_dense_true.T, c='gray', s=8, alpha=0.3, label='triangulated with true poses')\n", "ax.scatter(*recon_dense_estimated.T, c=point_colors, s=25, label='reconstructed (estimated pose)')\n", "\n", "# mark the two camera centers\n", "cam1_center = -R1.T @ t1\n", "cam2_center = -R_true.T @ t_true\n", "ax.scatter(*cam1_center, c='blue', s=80, marker='^', label='camera 1')\n", "ax.scatter(*cam2_center, c='green', s=80, marker='^', label='camera 2')\n", "\n", "ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')\n", "ax.legend(fontsize=8)\n", "ax.set_title('Structure from motion: recovered 3D points and camera poses')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d4820c99", "metadata": {}, "source": "### From sparse to dense\n\nEven with the loosened thresholds above, feature matching only yields a *sparse* cloud — one 3D point per distinctive image keypoint. As shown above, the sparse cloud is often not very satisfying. But with the camera poses now known, we're no longer limited to distinctive keypoints: **multi-view stereo (MVS)** estimates a depth for *every* pixel in *every* view, by sweeping a plane-hypothesis or patch through space (a generalization of Lesson 22's stereo matching from a rectified pair to arbitrarily-posed calibrated cameras). Those per-view depth maps are then fused — inconsistent estimates filtered out, the rest merged — into the dense colored point cloud (or mesh) that toolkits like COLMAP (Schönberger & Frahm, 2016) produce." }, { "cell_type": "markdown", "id": "a5fd9ef5", "metadata": {}, "source": [ "### From two views to many: incremental SfM and bundle adjustment\n", "\n", "Real reconstructions rarely stop at two views. The standard recipe sequentially incorporates new images in 3 steps: 1) match features between the new image and one or more existing images; 2) compute the new camera pose relative to the existing 3D point cloud; and 3) update the 3D point cloud with the new matches by triangulation. Step 2 requires solving Perspective-n-Point (i.e., given 2D-3D correspondences and $K$, recover that camera's pose via `cv2.solvePnP`).\n", "\n", "As views are accumulated, errors compound — pose corrupts the triangulated points, which then corrupt the next recovered pose, and so on. **Bundle adjustment** fixes this by refining all camera poses and 3D points *together* in one large nonlinear least-squares optimization that minimizes the total reprojection error summed across every point in every view at once. Bundle adjustment is arguably the single most important step separating a simple two-view pipeline like the one in this lesson from a real SfM system (COLMAP, etc.)." ] }, { "cell_type": "markdown", "id": "f5507346", "metadata": {}, "source": [ "### Exercises\n", "\n", "1. Add extra synthetic pixel noise (e.g. `rng.normal(0, 1.0, x1_in.shape)`) on top of the already-real detections in `x1_in`/`x2_in` before triangulating. How much does the reconstruction error grow, and does it grow uniformly, or worse for points farther from the cameras (Lesson 22's disparity-depth relationship: distant points produce smaller, noisier parallax)?\n", "2. Look up `cv2.solvePnP`'s signature and sketch (in words) how you'd extend this notebook to a third view: which 3D points from the `view1`/`view2` reconstruction would you match into the new image, and how would `solvePnP`'s inputs and outputs map onto the pieces we already have (`K`, matched 2D points, and the existing 3D point cloud)?\n", "3. This pair (views 0 and 2) has a deliberately generous baseline for reliable matching. Swap in views 0 and 1 instead (a much shorter baseline — you'll need their `.npz` camera parameters from the same BlendedMVS scene). How does the shorter baseline affect the number of RANSAC inliers, and the accuracy of the recovered pose and reconstruction, compared to the wider-baseline pair used above?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }