Using Python 3.9 or newer, create a virtual environment as follows:
python3 -m venv venv
source venv/bin/activate
Install this homework's dependencies:
pip install -r requirements.txt
You can now open the project directory in VS Code. Within VS Code, open the command palette (⌘ command ⇧ shift P), run Python: Select Interpreter, and choose the virtual environment you created in the previous steps.
For the best editing experience, install the following VS Code extensions:
- Python (
ms-python.python) - Pylance (
ms-python.vscode-pylance) - Black Formatter (
ms-python.black-formatter) - Ruff (
charliermarsh.ruff)
In this problem set, you'll be implementing a struction-from-motion pipeline (SfM) from the ground up, from camera calibration to fundamental matrix estimation. You'll be testing your pipeline either in a simulated environment or in the real world to create a 3D point cloud of an environment from multi-view images.
This problem set was created by Christian Arnold (cmarnold@mit.edu) and Isabella Yu (iyu@mit.edu).
Images that you will create will be save to their respective problem folder in outputs/p[x]_problem_name
The most important part of any SfM pipeline is calibrating the camera. This is commonly done with a chessboard object containing a specified number of squares. The best calibration comes from chessboards that space the entirety of the image plane, so chessboard images are not necessarly square like a normal chessboard. In this exercise, you will identify the edges, or the contours of the image. This is a critical step in calibrating the camera. There are many possible implementations for edge detectors, we only require that you not use an existing utility like opencv to detect contours, or other similiar libraries.
Complete the ContourImage() class p1_edge_identification.py to prepare the image for processing. Specifically, complete definitions for binarize(self) and fill_border(self). This is a typical step in edge detection in imagery, so if you are confused, we encourage you to review online resources.
Complete the find_contours(binary_image, foreground=1) function definition in p1_edge_identification.py, then run python -m src.p1_edge_identification. The output of this function should be a one-dimensional list of (u, v) coordinates of every edge in the image. Your final output will look similar to the image below.
From the edge contours of the chessboard image, it is trivial, yet time-consuming, to then form quadrilaterals and detect the corners of each square. We have instead used the existing OpenCV implementation of cv2.findChessboardCorners. With the chessboard corners, we can now recover the intrinsic properties of the camera, and correct for any distortions. There are many methods of recovering camera calibration, and each involves assumptions about the camera, as such, no implementation will be exactly the same. We provide a wrapper for the OpenCV implementation of cv2.calibrateCamera.
In previous psets, you have already been introduced to the Intrinsic Camera Matrix, K,
Where f_x, f_y are the focal ratios for that axis and c_x, c_y are the u, v coordinates for the center of the image.
Imagine a pinhole camera with a small aperture (“pinhole”) at the origin of our coordinate system. The camera is oriented so that the optical (principal) axis is the Z-axis, pointing forward. The image plane is placed at Z = f (focal length f ) in front of the pinhole. A pinhole camera creates a perfect image with no distortions where f_x = f_y = f_r * min(image_height, image_width) and f_r = 1 / 2tan(FoV/2) where the FoV (Field of View) is in radians. Assume the FoV is 45 degrees. Calculate an intrinsic matrix for the image found in data/p1_edge_identification/chessboard.png and assign the value of ideal_instrinsic_matrix in p2_calibrate_camera.py.
In p2_calibrate_camera.py, define the interior dimensions (number of squares - 1) of chessboard_size. You should see an image similar to the one below with the corners drawn on the chessboard.
To calibrate the camera, we need to compare the image corners with the corresponding points of a known 3D object. In this case, it is the chessboard. We can assume the surface of the chessboard is perfectly flat, so all z-values of the object points should be 0. Complete the get_3D_object_points(chessboard_size) in p2_calibrate_camera.py which defines the corresponding 3D points for each corner. You can assume the first point is (0, 0, 0) and the last point is (chessboard_size[0]-1, chessboard_size[1]-1, 0). The 3D points should be in the same order as the recovered chessboard corners.
Once you have correctly defined the 3D points, you should be able to recover the intrinsic matrix and estimated distortion coefficients. We have artificially distorted the image, so the recovered intrinsic matrix may be different from the ideal matrix you computed in 2.a.
Note: All methods of camera calibration are providing an over-fitted guess, thus the recovered distorition will not be the same as the actual distortion applied to the image.
The most common form of image distortion is radial distortion. Every real camera has some amount of radial distortion that must be corrected for Multi-View Geometry.
Distortion Coefficients recovered by OpenCV are defined as
The below equation applies a distortion.
where
Your task is to undo the distortion of the original image using the recovered intrinsic matrix and distortion coefficients by completing the definition of undistort_image(image, camera_matrix, dist_coeffs) in p2_calibrate_camera.py.
Run python -m src.p2_camera_calibration to recover the camera parameters from simulation.
Do not use OpenCV functions. If you are stuck, you may use cv2.getOptimalNewCameraMatrix and cv2.undistort for half credit. If you do use OpenCV, it is helpful to break the problem down into two parts:
- Recover the Optimal Camera Matrix, and
- Undistort the image.
In this part, you will work on the geometric relationship between two views of a scene. Your goal is to compute the fundamental matrix using two different methods and then use the matrix to draw epipolar lines that explain how points in one image constrain the location of corresponding points in the other image.
Implement the lstsq_eight_point_alg in p3_fundamental_matrix.py. This method forms a system of equations based on corresponding points in two images. By stacking the constraints into a matrix (denoted as W), the method uses singular value decomposition (SVD) to solve for the fundamental matrix F that satisfies:
After obtaining an initial solution, enforce that F has rank 2 by zeroing out the smallest singular value. Do not use OpenCV functions in this section.
Implement the normalized_eight_point_alg in p3_fundamental_matrix.py. In the normalized eight-point algorithm, you first normalize the coordinates of the input points by translating the points so that their centroid is at the origin and scaling the points so that the mean squared distance from the origin is 2.
The normalization is done with transformation matrices
The normalization further refines
Now implement compute_epipolar_lines. Once you have a fundamental matrix, you can compute the corresponding epipolar line for any given point in one image. The epipolar line is computed by multiplying the fundamental matrix with the homogeneous coordinate of a point, where l = F.dot(p).
Run python -m src.p3_image_rectification to produce the resulting images. You should now see that all epipolar lines pass through or close to the defined points.
In this part, you will work on aligning (rectifying) two images so that their epipolar lines become horizontal. Rectification is a critical step for stereo vision, as it simplifies correspondence search between images.
The epipoles are the points of intersection of the camera's optical axis with the image plane. They serve as the "center of projection" for the corresponding epipolar lines in each image.
Implement compute_epipole in p4_image_rectification.py. Compute the epipolar lines by multiplying the fundamental matrix with the points from the other image. Use singular value decomposition (SVD) to solve for the epipole that satisfies the homogeneous equation l.dot(e) = 0, where l is the epipolar line and e is the epipole. Normalize the epipole so that the last coordinate is 1. Do not use OpenCV functions in this section.
Once you have the epipoles, the next step is to compute the homographies that will align (rectify) the images.
Implement the function compute_matching_homographies.
For the first image, use the computed
For the second image, translate the image so that its center is at the origin. Rotate the image so that the epipole is aligned with the x-axis. Apply a projective transformation that sends the epipole to infinity. Combine these steps into a single homography
Solve for the affine component using a least-squares fit (np.linalg.lstsq) that relates the transformed matching points in the two images. Do not use OpenCV functions in this section.
After computing the homographies, use them to transform the images into a rectified coordinate system where the epipolar lines are horizontal.
Implement the function compute_rectified_image(im, H) which applies the homography H to an image. First, compute the new coordinates for each pixel using the homography. Determine the new image bounds and calculate the offsets. Use the inverse homography to map every point from the rectified image back to the original image to sample the pixel values. Return the rectified image and the offsets, which will be used to align the two images.
Once rectified, corresponding points in the two images should lie on the same horizontal line, which can be verified by the output image. Do not use OpenCV functions in this section.
Now, run python -m src.p4_image_rectification get the rectified images of our test object. The rectification output should look like this:
There have been many algorithms developed for feature detection over the years, with each iteration improving on the last. SIFT, SURF, and ORB are all popular implementations of feature detectors that common points between different images using keypoints and descriptors for those points. ORB is an improvement of SIFT and SURF that tends to produce the most matched features.
Implement find_matches using cv2.SIFT_create() to identify features between two images. Then use Fast Linear Nearest Neighbors (FLANN) matching from cv2.FLANNMatcher with FLANN_INDEX_KDTREE = 1, index_params = dict(algorithm = FLANN_INDEX_KDTREE, trees = 5), and search_params = dict(checks=50) to match the possible features.
Filter for good_matches by using a ratio threshold of 0.75.
This is the last component of a basic Structure from Motion pipeline. More complex versions of SfM have further refining steps like dense reconstruction and texturing, but we will stop at a sparse point cloud. We'll be constructing a point cloud of a model of the Arc de Triomphe, which is found in data/p5_3D_reconstruction/arc_de_triomphe
Implement recover_fundamental_matrix
Leverage OpenCV’s RANSAC-based method to compute the fundamental matrix from matched feature points between two images. Use a threshold of 1.0 and confidence of 0.99.
Implement compute_essential_matrix
Convert the fundamental matrix into the essential matrix using the intrinsic camera matrix. The essential matrix relates corresponding normalized image points and is used to derive the relative camera pose.
Implement estimate_initial_RT.
Use singular value decomposition (SVD) of the essential matrix to obtain two candidate rotations and two candidate translations. These are potential relative poses between the camera views.
Implement find_best_RT.
To resolve ambiguities, triangulate the 3D points for each candidate pose. The correct pose is the one that yields the maximum number of points with positive depth in both camera coordinate systems. You are encouraged to use cv2.triangulatePoints to solve for which combination of candidate rotations and translations results in a positive depth.
Implement get_identity_projection_matrix and get_relative_projection_matrix. The identity projection should be the dot product of the camera matrix and [I|0]. P2 should similarly be the the dot product of the camera matrix and [R|T].
With the projection matrixes defined, we can now triangulate the points relative to the projections. This yeilds a 4D homogenous point, which we convert to 3D points. You can see the 3D points in the matplotlib plot.
The singular values of the fundamental matrix are also listed as a sanity check. F should have a rank of 2. If the second value in value is very close to zero, the solution may be degenerate, which would result in a linear point cloud.
Now, you can run python -m src.p5_3D_reconstruction to get a point cloud from stereo views of the Arc de Triomphe. Run with the --setup flag if you want to generate new images of the Arc de Triomphe (may take a while!).
Expected SIFT output:
Expected point cloud output:
Now, we can use all of the previous functions to get a point cloud from more than two images. This might require bundle adjustment for jointly optimizing camera parameters and the point cloud, as well as defining some correspondences manually, to get a reasonable looking point cloud. Implement bundle adjustment, reprojection_error, and main in p6_SFM_pipeline.py to get a point cloud of the statue dataset env.p5.statue_images from 5 or more views. We will only be grading your implementation of bundle_adjustment and reprojection_error. Then run python -m src.p6_SfM_pipeline to produce a point cloud.
This was our output point cloud from the statue dataset (using some handcrafted points):
Feel free to run camera calibration on your own phone using a printed-out checkerboard if you want to capture real-world images! We recommend printing out the checkerboard intead of taking pictures of the checkerboard on your screen to reduce moire effects. The checkerboard should cover the entire FoV of the camera.
Parts of this pset were inspired by CS231A, taught by Jeannette Bohg and Silvio Savarese in spring 2024. The Arc de Triomphe model can be found here: https://rigmodels.com/model.php?view=Arc_de_Triomphe_3d_model__65937fd27de647c0a8ac99ce8275c03e. It has a royalty-free license.