Image source: BlendedMVS (CC BY 4.0)
", "metadata": {} }, { "cell_type": "markdown", "id": "20a17371", "source": "## Step 1: feature matching\n\nSIFT + ratio test (Lesson 19), run directly on the two photos.", "metadata": {} }, { "cell_type": "code", "id": "323b120e", "source": "gray1 = cv2.cvtColor(view1, cv2.COLOR_BGR2GRAY)\ngray2 = cv2.cvtColor(view2, cv2.COLOR_BGR2GRAY)\n\nsift = cv2.SIFT_create()\nkp1, des1 = sift.detectAndCompute(gray1, None)\nkp2, des2 = sift.detectAndCompute(gray2, None)\n\nbf = cv2.BFMatcher()\nraw_matches = bf.knnMatch(des1, des2, k=2)\ngood_matches = [m for m, n in raw_matches if m.distance < 0.75 * n.distance]\n\nx1 = np.float32([kp1[m.queryIdx].pt for m in good_matches])\nx2 = np.float32([kp2[m.trainIdx].pt for m in good_matches])\n\nprint(f'keypoints: {len(kp1)} (view 1), {len(kp2)} (view 2)')\nprint(f'good matches after ratio test: {len(good_matches)}')\n\nmatch_vis = cv2.drawMatches(view1, kp1, view2, kp2, good_matches, None,\n flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)\nplt.figure(figsize=(10, 4))\nplt.imshow(cv2.cvtColor(match_vis, cv2.COLOR_BGR2RGB))\nplt.title('Feature matches')\nplt.axis('off')\nplt.show()", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "313721d7", "metadata": {}, "source": "## Step 2: robust essential matrix and pose recovery\n\nFeed the matched points through the same robust pipeline built in Lessons 22 and 26 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 26, `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 23) and the fundamental matrix (Lesson 26):\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$ (smallest-singular-value right singular vector) 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\ninliers = inlier_mask.ravel().astype(bool)\nx1_in, x2_in = x1[inliers], x2[inliers]\n\nrecon_true_cameras = triangulate_dlt(P1, P2_true, x1_in, x2_in)\n\ncv_points_4d = cv2.triangulatePoints(P1, P2_true, x1_in.T, x2_in.T)\ncv_points_3d = (cv_points_4d[:3] / cv_points_4d[3]).T\nprint(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\ndepth1 = recon_true_cameras[:, 2]\ndepth2 = (R_true @ recon_true_cameras.T + t_true.reshape(3, 1))[2]\nvalid = (depth1 > 0) & (depth2 > 0)\nx1_in, x2_in, recon_true_cameras = x1_in[valid], x2_in[valid], recon_true_cameras[valid]\nprint(f'kept {valid.sum()} / {len(valid)} inliers after cheirality check')\n\ndef 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\nerr1 = reprojection_error(P1, recon_true_cameras, x1_in)\nerr2 = reprojection_error(P2_true, recon_true_cameras, x2_in)\nall_err = np.concatenate([err1, err2])\nprint(f'reprojection error using TRUE camera poses: mean {all_err.mean():.3f} px, max {all_err.max():.3f} px')" }, { "cell_type": "markdown", "id": "3d4fc872", "metadata": {}, "source": "## Putting it together: reconstruction using *estimated* pose\n\nNow 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 (the one piece of information triangulation from image correspondences alone can never supply — here read off from BlendedMVS's ground truth, standing in for whatever real-world scale reference you'd use in practice) 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)])\nreconstruction_unit_scale = triangulate_dlt(P1, P2_estimated, x1_in, x2_in)\n\ntrue_baseline = np.linalg.norm(t_true)\nreconstruction_metric = reconstruction_unit_scale * true_baseline\n\nerror = np.linalg.norm(reconstruction_metric - recon_true_cameras, axis=1)\nprint(f'true baseline: {true_baseline:.3f}')\nprint(f'mean 3D reconstruction error (vs. triangulation from true poses): {error.mean():.3f}')\nprint(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 21) instead of two arbitrary independent cameras." }, { "cell_type": "markdown", "id": "df0ab549", "metadata": {}, "source": "### Visualizing the reconstruction\n\nFor 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", "id": "42b6c294", "source": "sift_dense = cv2.SIFT_create(contrastThreshold=0.01)\nkp1d, des1d = sift_dense.detectAndCompute(gray1, None)\nkp2d, des2d = sift_dense.detectAndCompute(gray2, None)\nraw_dense = bf.knnMatch(des1d, des2d, k=2)\ngood_dense = [m for m, n in raw_dense if m.distance < 0.85 * n.distance]\n\nx1d = np.float32([kp1d[m.queryIdx].pt for m in good_dense])\nx2d = 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)\ndense_inliers = dense_mask.ravel().astype(bool)\nx1d, x2d = x1d[dense_inliers], x2d[dense_inliers]\n\nrecon_dense_true = triangulate_dlt(P1, P2_true, x1d, x2d)\nrecon_dense_estimated = triangulate_dlt(P1, P2_estimated, x1d, x2d) * true_baseline\n\n# cheirality check again, same as for the sparse set above\ndepth1 = recon_dense_true[:, 2]\ndepth2 = (R_true @ recon_dense_true.T + t_true.reshape(3, 1))[2]\nvalid = (depth1 > 0) & (depth2 > 0)\nx1d, recon_dense_true, recon_dense_estimated = x1d[valid], recon_dense_true[valid], recon_dense_estimated[valid]\n\nrgb1 = cv2.cvtColor(view1, cv2.COLOR_BGR2RGB)\npx = np.clip(x1d[:, 0].astype(int), 0, rgb1.shape[1] - 1)\npy = np.clip(x1d[:, 1].astype(int), 0, rgb1.shape[0] - 1)\npoint_colors = rgb1[py, px] / 255.0\n\nprint(f'{len(x1d)} points in the denser cloud (vs. {len(x1_in)} used for the numeric check above)')", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "execution_count": null, "id": "e9360d20", "metadata": {}, "outputs": [], "source": "fig = plt.figure(figsize=(7, 6))\nax = fig.add_subplot(111, projection='3d')\nax.scatter(*recon_dense_true.T, c='gray', s=8, alpha=0.3, label='triangulated with true poses')\nax.scatter(*recon_dense_estimated.T, c=point_colors, s=25, label='reconstructed (estimated pose)')\n\n# mark the two camera centers\ncam1_center = -R1.T @ t1\ncam2_center = -R_true.T @ t_true\nax.scatter(*cam1_center, c='blue', s=80, marker='^', label='camera 1')\nax.scatter(*cam2_center, c='green', s=80, marker='^', label='camera 2')\n\nax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')\nax.legend(fontsize=8)\nax.set_title('Structure from motion: recovered 3D points and camera poses')\nplt.show()" }, { "cell_type": "markdown", "id": "d4820c99", "source": "### From sparse to dense\n\nEven loosened, feature matching only ever gives a *sparse* cloud — one point per distinctive keypoint, with large gaps over smooth, low-texture regions. Once SfM has recovered the camera poses, though, every pixel becomes fair game: **multi-view stereo (MVS)** sweeps a plane-hypothesis or patch (Lesson 21's stereo matching, generalized from a rectified pair to arbitrarily-posed calibrated cameras) through space to estimate a depth for *every* pixel in *every* view, then fuses those per-view depth maps — filtering out inconsistent estimates and merging the rest — into the dense colored point cloud (or mesh) that toolkits like COLMAP or, fittingly, BlendedMVS's own reconstruction pipeline produce.", "metadata": {} }, { "cell_type": "markdown", "id": "f5507346", "metadata": {}, "source": "### Exercise\n\n1. 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 21's disparity-depth relationship: distant points produce smaller, noisier parallax)?\n2. This lesson only used 2 of the scene's 131 available views. Real structure-from-motion pipelines (e.g. COLMAP) add many more views incrementally: each new image is *resected* against the already-reconstructed 3D points (a problem called Perspective-n-Point, or PnP — given 2D-3D correspondences and $K$, solve for that camera's pose) via `cv2.solvePnP`, then its new points are triangulated against the growing reconstruction. Look up `cv2.solvePnP`'s signature and sketch (in words) how you'd extend this notebook to a third view.\n3. 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 }