From 2D Images to 3D Worlds: The Linear Algebra of Structure from Motion

Bach Pham

Abstract

This paper explores the essence of Structure from Motion (SfM), a technique in computer vision that can reconstruct 3D scenes from some 2D images with good overlap. The paper aims to apply linear algebra, especially Singular Value Decomposition (SVD), to deeply explore the mathematics behind SfM. I focus mostly on the Fundamental Matrix and the Eight-Point Algorithm. Specifically, I demonstrate how geometric constraints are transformed into a homogeneous system of linear equations ($Ax = 0$) and prove why the Singular Value Decomposition provides the optimal solution to this problem.

Introduction: The Inverse Problem

Sometimes we only have a handful of images of a scene, but we really want to know what the whole scene looks like. People have thought of combining images to create panoramas using feature matching, but this is still just a larger image. We want to move from a flat 2D image to a 3D structure.

In the 1970s, structure from motion was developed. It has been the ground to many modern techniques to reconstruct 3d scene, especially all kinds of Neural Radiance Field methods (NeRFs) and Gaussian Splatting.

In rendering (forward problem), we project 3D coordinates onto a 2D plane, which causes the loss of depth information. On the other hand, Structure from Motion (SfM) tries to reverse this by using the parallax effect which is the relative motion of points between different viewpoints.

A prominent implementation of this is Google Map. It has been using SfM all the time until 2023 when Google released the "immersive view" option on Google Maps which uses another technology built on top of SfM: NeRF. The NeRF paper was released in 2020 by Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng. Later on in 2024, Google released a research paper about a hybrid technology (combining the high quality of NeRF and the ultra fast rendering of Gaussian Splatting): RadSplat

The SfM pipeline generally consists of:

1. Feature Matching (SIFT) and Outlier Rejection (RANSAC).
2. Estimating the Fundamental Matrix ($F$) (The focus of this paper).
3. Estimating the Essential Matrix ($E$).
4. Triangulation and Bundle Adjustment.

Epipolar Geometry

Before we get to the Fundamental Matrix, we have to know what Epipolar geometry is.

The Coplanarity Constraint

Consider a point $\bX$ in 3D space viewed by two cameras centered at $\bC$ and $\bC'$.

• Let $\bx$ be the projection of $\bX$ onto the first image.
• Let $\bx'$ be the projection of $\bX$ onto the second image.

As illustrated in Figure 1, the points $\bC$, $\bC'$, and $\bX$ form a triangle. Consequently, the vectors $\overrightarrow{\bC\bx}$, $\overrightarrow{\bC\bC'}$, and $\overrightarrow{\bC'\bx'}$ must lie on the same plane, called the epipolar plane. And since they are on the same plane, their scalar triple product is zero:

Epipoles and Epipolar Lines

If we know a point $\bx$ in the above image, we cannot determine its exact depth. However, the ray from $\bC$ through $\bx$ looks like a line when viewed from the second camera. This is the Epipolar Line ($\mathbf{l}'$). The matching point $\bx'$ must lie somewhere on this line. This constraint significantly reduces the search for matches because we just need to search on a line and do not need to search for the entire image anymore.

The Fundamental Matrix $F$

The Fundamental Matrix $F$ is the algebraic representation of epipolar geometry. It is a $3 \times 3$ matrix. The relationship between a point in image 1, $\mathbf{x}_i = [x_i, y_i, 1]^T$, and a point in image 2, $\mathbf{x}'_i = [x'_i, y'_i, 1]^T$, is given by the Longuet-Higgins equation:

Deriving the System of Equations ($Ax=0$)

To find $F$, we need to solve for its 9 entries. We can write the matrix product explicitly for a single pair of points:

Expanding this multiplication results in a scalar equation. Grouping the terms to isolate the coefficients of $F$, we obtain:

If we have $m$ correspondences (where $m \geq 8$), we can stack these linear equations into a large matrix system. By defining a vector $\mathbf{f}$ that contains the entries of the matrix $F$, we can rewrite the system as:

This can be written compactly as $A\mathbf{f} = 0$.

Since real-world data has noise, $A\mathbf{f}$ will never be exactly zero. We instead try to find the $\mathbf{f}$ that makes the result as close to zero as possible. We also fix the length of $\mathbf{f}$ to 1 so the solution isn't just a vector of all zeros.

Mathematical Proof: Solving via SVD

Understand SVD geometrically

We have already known that $A = U \Sigma V^T$. Geometrically, multiplying a vector $\mathbf{x}$ by $A$ happens in three steps:

1. Rotate ($V^T$): The matrix $V^T$ is orthogonal. It takes the input vector $\mathbf{x}$ and rotates it to a new coordinate system.
2. Stretch ($\Sigma$): The matrix $\Sigma$ is diagonal. It stretches the vector along the coordinate axes. The amount of stretching is determined by the singular values $\sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_n \geq 0$.
3. Rotate ($U$): The matrix $U$ is orthogonal. It rotates the stretched vector one final time.

We want the final result $A\mathbf{x}$ to be as small as possible.

Proof

Step 1: Ignore the Final Rotation ($U$) We have:

The matrix $U$ is orthogonal. A property of orthogonal matrices is that they act like rotation or reflection matrix. That means they change the direction of a vector but not the length. Therefore, for any vector $\mathbf{z}$, $\| U\mathbf{z} \| = \| \mathbf{z} \|$.

Step 2: Change Variables. Let's define a temporary vector $\mathbf{y} = V^T \mathbf{x}$. Since $V$ is orthogonal, $V^T$ is also orthogonal. This means $V^T$ just rotates $\mathbf{x}$. Since we are constrained to unit vectors ($\| \mathbf{x} \| = 1$), and rotation preserves length, our new vector $\mathbf{y}$ is also a unit vector ($\| \mathbf{y} \| = 1$).

Our problem is now much simpler: Minimize $\| \Sigma \mathbf{y} \|$.

Step 3: Analyze the Stretching ($\Sigma$)

The squared length of this resulting vector is the sum of its squared components:

We want to make this sum as small as possible. We know the singular values are sorted in descending order: $\sigma_1$ is the largest (maximum stretch) and $\sigma_n$ is the smallest (minimum stretch). We also know that $y_1^2 + y_2^2 + \dots + y_n^2 = 1$ (because $\mathbf{y}$ is a unit vector).

To get the smallest total sum, we should put all of our weight into the term multiplied by the smallest number.

• If we choose $y_1=1$, the result is $\sigma_1^2$ (The Maximum).
• If we choose $y_n=1$, the result is $\sigma_n^2$ (The Minimum).

Therefore, the optimal vector $\mathbf{y}$ is all zeros, except for a 1 in the last position:

Step 4: Convert back to $\mathbf{x}$. We found the best $\mathbf{y}$, but we need to find the original $\mathbf{x}$. We defined $\mathbf{y} = V^T \mathbf{x}$. Therefore, x is

By the rules of matrix-vector multiplication, this operation selects the last column of the matrix $V$.

Conclusion: The vector $\mathbf{x}$ that minimizes the error is the last column of $V$.

Rank Enforcement

After using the theorem above, we have a solution for $\mathbf{f}$ and reshape it into a $3 \times 3$ matrix $F$. However, we have one problem.

The Rank 2 Constraint

In our geometric setup, the Fundamental Matrix maps a point to a line. For this to work mathematically, the determinant of $F$ must be zero. In Linear Algebra terms, $F$ must be Rank 2.

However, because our input images have noise (pixels aren't perfect), the $F$ we calculated in Section 4 will likely be Rank 3 (Full Rank). It will satisfy the algebraic least squares, but it won't satisfy the strict geometry (the epipolar lines won't intersect perfectly at one point).

The Eckart-Young-Mirsky Theorem

We need to find a new matrix $F'$ that is Rank 2 but is as close as possible to our calculated $F$. The Eckart-Young-Mirsky theorem tells us exactly how to do this:

1. Decompose: Take the SVD of our noisy matrix $F$: $F = U \Sigma V^T$.
2. Inspect $\Sigma$: The diagonal matrix will look like this:
   $$\Sigma = \operatorname{diag}(\sigma_1, \sigma_2, \sigma_3)$$
   Because of noise, $\sigma_3$ will be a small non-zero number (e.g., 0.005). This tiny number is what makes the matrix Rank 3.
3. Truncate: We manually set the smallest singular value to zero.
   $$\Sigma' = \operatorname{diag}(\sigma_1, \sigma_2, 0)$$
4. Reconstruct: We multiply the matrices back together: $F' = U \Sigma' V^T$.

This result, $F'$, is the Fundamental Matrix used in Structure from Motion. Of course this is just the very basic idea. We then need to use RANSAC again for outlier rejection to find best F, then estimate the essential matrix, then the triangulation, then PnP, and a lot of optimizations. However, I will not include those in this paper because we only focus on the core math, not the programming aspect.

Conclusion

SfM is a very interesting application of Linear Algebra. The math behind is so beautiful, especially for the minimization problem. It is also very convenient when we can just decompose, force an entry in the diagonal matrix to be zero, and multiply everything back. After taking linear algebra course, I finally understand deeply what is happening behind the libraries I imported in my code.

References

[1] University of Maryland. (2022). CMSC733 Project 3: Build your own SfM! Retrieved from https://cmsc733.github.io/2022/proj/p3/