HW3: 3D reconstruction

python q1.py -a True

python q1.py -b True

python q2.py -p 7 -o hydrant

python q3.py

Q1: 8-point and 7-point algorithm (40 points)

(A1) F matrix using 8-point algorithm (15 points)

1. Brief explanation of your implementation.

1. Data Preprocessing

Normalize the point coordinates to improve numerical stability. The normalization steps are as follows:

• Calculate the centroid and average distance of each set of points, and use these results to construct the normalization matrix T.
• Transform the original point coordinates into a new coordinate space using the normalization matrix T.

2. Construct Matrix A

In the normalized coordinate space, construct matrix A using each pair of matched points in the form:

    A = [x'1 x1   x'1 y1   x'1   y'1 x1   y'1 y1   y'1   x1   y1   1]

where (x₁, y₁) and (x'₁, y'₁) are the coordinates of corresponding points.

3. Solve for the Fundamental Matrix F

Perform Singular Value Decomposition (SVD) on matrix A and select the eigenvector corresponding to the smallest singular value (the null space of A) to reconstruct F:

• F is the null space vector of A, reshaped into a 3 × 3 matrix.
• Enforce the rank-2 constraint on F (a property of the fundamental matrix) by setting the smallest singular value to zero after SVD.

4. Denormalization

The obtained F is in the normalized coordinate space, so it is transformed back to the original coordinate system as follows:

F = T2T F T1

where T₁ and T₂ are the respective normalization matrices.

5. Compute Epipolar Lines

Using the fundamental matrix F and points in one image, calculate the corresponding epipolar lines in the other image. The epipolar line satisfies:

l' = F ⋅ p

where p is a point in image 1, and l' is the epipolar line in image 2.

Name	F-matrix visualizations
Bench
Bench
remote
remote

Fundamental Matrices

• bench
  • \(\begin{bmatrix}-1.19265833e-07 & -2.53255522e-06 & 2.07584109e-04 \\ -3.96465204e-07 & 2.27096267e-07 & 2.48867750e-02 \\ -1.06890748e-03 & -2.30052117e-02 & 1.00000000e+00 \end{bmatrix}\)
• remote
  • \(\begin{bmatrix} 7.19006698e-07 & -1.90583125e-07 & 2.36215166e-03 \\ -1.88159055e-06 & -8.60510729e-07 & -8.45685523e-03 \\ -3.98058321e-03 & 9.46500248e-03 & 1.00000000e+00 \end{bmatrix}\)

(A2) E matrix using 8-point algorithm (5 points)

1. Brief explanation of your implementation.

Compute the Essential Matrix E

1. Calculate the Essential Matrix E

The essential matrix E is calculated using the fundamental matrix F and the intrinsic matrices K₁ and K₂ of the cameras. The formula is:

E = K2T ⋅ F ⋅ K1

where:

• K₁ and K₂ are the intrinsic matrices of the two cameras.
• F is the fundamental matrix.

2. Normalization

To ensure consistency in the scale of the essential matrix, normalize E so that the last element of E equals 1:

E = E / E[2, 2]

2. Provide your estimated E.

• bench
  • \(\begin{bmatrix} -0.41524274 & -8.81748894 & -2.15055808 \\ -1.3803559 & 0.79067134 & 69.08265036 \\ -3.58890876 & -68.47438701 & 1.00000000 \end{bmatrix}\)
• remote
  • \(\begin{bmatrix} 2.5298064 & -0.67056178 & 7.33094837 \\ -6.62032749 & -3.02768466 & -28.29861665 \\ -14.50071092 & 26.41824781 & 1.00000000 \end{bmatrix}\)

(B) 7-point algorithm (20 points)

1. Brief explanation of your implementation.

1. Data Preprocessing

Normalize the point coordinates to improve numerical stability:

• For each set of corresponding points, calculate the centroid and average distance, and use these values to construct a normalization matrix T.
• Transform the original point coordinates into a new coordinate space using the normalization matrix T.

2. Construct Matrix A

In the normalized coordinate space, construct matrix A using each pair of matched points in the form:

    A = [x'1 x1   x'1 y1   x'1   y'1 x1   y'1 y1   y'1   x1   y1   1]

where (x₁, y₁) and (x'₁, y'₁) are the coordinates of corresponding points.

3. Solve for Candidate Fundamental Matrices F

Perform Singular Value Decomposition (SVD) to obtain the null space of matrix A:

• Use SVD to get the last two singular vectors of A, f₁ and f₂, each reshaped into a 3 × 3 matrix, representing two possible solutions for the fundamental matrix.
• Use the linear combination F = α f₁ + (1 - α) f₂ and substitute it into the constraint det(F) = 0 to form a cubic polynomial in α.

4. Solve the Polynomial Equation

Solve the cubic polynomial to obtain possible real roots for α, with each real root providing a candidate fundamental matrix F.

5. Denormalization

Transform each candidate fundamental matrix from the normalized coordinate space back to the original coordinate system:

F = T2T F T1

where T₁ and T₂ are the normalization matrices.

Name	F-matrix visualizations
hydrant
hydrant
ball
ball

Fundamental Matrices

• hydrant
  • \(\begin{bmatrix} -3.35078059e-09 & -2.91545005e-06 & -6.74891733e-03 \\ 3.51682683e-06 & -8.84237457e-07 & -1.50044163e-02 \\ 9.09398538e-03 & 1.48748649e-02 & 1.00000000e+00 \end{bmatrix}\)
• ball
  • \(\begin{bmatrix} 1.06836663e-07 & -3.30854241e-06 & -1.37397509e-03 \\ 4.79268901e-06 & 1.63305018e-07 & -1.67456794e-02 \\ 1.12974017e-03 & 1.64136613e-02 & 1.00000000e+00 \end{bmatrix}\)

Q2: RANSAC with 7-point and 8-point algorithm (20 points)

1. Brief explanation of your RANSAC implementation and criteria for considering inliers.

1. Overview of the RANSAC Algorithm

RANSAC (Random Sample Consensus) is an iterative method used to find inliers and estimate model parameters from data that contains noise. For fundamental matrix estimation, the goal is to obtain a robust estimation of F by removing outliers.

2. Random Sampling and Model Fitting

In each iteration:

• Randomly select 7 or 8 points from the set of matched points (depending on whether the 7-point or 8-point algorithm is used).
• Calculate the fundamental matrix F using the selected points through the 7-point or 8-point algorithm. In your code, the compute_F function selects the corresponding algorithm based on the number of points.

3. Error Calculation and Inlier Counting

For the calculated fundamental matrix F:

• Calculate the epipolar error for all matched points (i.e., the epipolar line error), achieved through the get_epi_error function.
• Points with errors below a dynamic threshold (dynamic_threshold) are considered inliers, and the inlier count is recorded. If the current inlier count exceeds the previous best record, the best fundamental matrix and inlier count are updated.

4. Dynamic Threshold and Early Stopping

To improve stability and efficiency, the threshold is reduced every 10% of the total iterations, tightening the inlier criteria.

If the inlier count remains stable for multiple rounds (up to early_stop_rounds), the iteration process terminates early.

5. Save the Best Fundamental Matrix F

The fundamental matrix F with the highest inlier count is returned as the best estimate.

Name	Algorithm	F-matrix visualizations
bench	8-Point
bench	8-Point
remote	8-Point
remote	8-Point
ball	8-Point
ball	8-Point
hydrant	8-Point
hydrant	8-Point
bench	7-Point
bench	7-Point
remote	7-Point
remote	7-Point
ball	7-Point
ball	7-Point
hydrant	7-Point
hydrant	7-Point

Fundamental Matrices

• 8 Point
  • bench
    • \(\begin{bmatrix} -6.53243663e-08 & -1.63821831e-06 & 1.53389833e-05 \\ -1.11114957e-06 & 2.82921707e-07 & 2.39057182e-02 \\ -9.19773080e-04 & -2.21329008e-02 & 1.00000000e+00 \end{bmatrix}\)
  • remote
    • \(\begin{bmatrix} 1.75660365e-06 & -4.36012735e-06 & 3.10230393e-03 \\ 4.27254212e-07 & 5.96576035e-07 & -8.92164955e-03 \\ -4.36166213e-03 & 9.33581306e-03 & 1.00000000e+00 \end{bmatrix}\)
  • ball
    • \(\begin{bmatrix} 3.57996410e-06 & 4.82772954e-06 & -3.60177320e-02 \\ -4.07477449e-06 & -3.50309464e-06 & -5.11016184e-02 \\ 4.76227305e-02 & 5.38368956e-02 & 1.00000000e+00 \end{bmatrix}\)
  • hydrant
    • \(\begin{bmatrix} 8.62317916e-09 & -8.22857408e-06 & 1.50256954e-04 \\ 1.01068463e-05 & -5.68964229e-07 & -2.03613343e-02 \\ -7.06139192e-04 & 2.07827792e-02 & 1.00000000e+00 \end{bmatrix}\)
• 7 Point
  • bench
    • \(\begin{bmatrix} -6.47824143e-07 & -1.15717448e-05 & 1.90263421e-03 \\ 7.48924387e-06 & -4.25597609e-06 & 3.19222235e-02 \\ -2.41042568e-03 & -2.74927906e-02 & 1.00000000e+00 \end{bmatrix}\)
  • remote
    • \(\begin{bmatrix} 1.68547776e-06 & -8.51910282e-06 & 3.86593182e-03 \\ 9.79177887e-06 & 5.49404065e-07 & -6.66720018e-03 \\ -7.02464122e-03 & 6.10542915e-03 & 1.00000000e+00 \end{bmatrix}\)
  • ball
    • \(\begin{bmatrix} -1.48290568e-07 & -3.46204352e-06 & -5.06111476e-03 \\ 4.07635331e-06 & -9.58646918e-07 & -1.35833710e-02 \\ 6.96238474e-03 & 1.35962310e-02 & 1.00000000e+00 \end{bmatrix}\)
  • hydrant
    • \(\begin{bmatrix} 3.57588730e-07 & 2.06762565e-06 & -4.04883975e-03 \\ -1.03328044e-06 & 8.77406599e-07 & -9.33808102e-03 \\ 3.93610666e-03 & 7.91780086e-03 & 1.00000000e+00 \end{bmatrix}\)

Name	Algorithm	% of inliers vs. # of RANSAC iterations
bench	8-Point
remote	8-Point
ball	8-Point
hydrant	8-Point
bench	7-Point
remote/td>	7-Point
ball	7-Point
hydrant	7-Point

Q3: Triangulation (20 points)

1. Brief explanation of your implementation.

1. Skew-Symmetric Matrix Construction

Convert each 2D point in the image to a skew-symmetric matrix. Given a point x = (x, y, 1)^T, the skew-symmetric matrix x_skew is:

    xskew = 
    [  0  -1   y ]
    [  1   0  -x ]
    [ -y   x   0 ]

This matrix allows us to represent the cross product as a matrix multiplication, facilitating constraint construction.

2. Constraint Matrix Calculation

For each point correspondence, calculate constraint matrices using the skew-symmetric representation:

    Constraint1 = x1, skew ⋅ P1
    Constraint2 = x2, skew ⋅ P2

where P₁ and P₂ are the camera projection matrices, and x₁ and x₂ are the corresponding points in the two images. These constraint matrices will be used to solve for the 3D point.

3. Triangulating 3D Points

For each pair of constraints Constraint₁ and Constraint₂:

• Stack the first two rows of each constraint matrix to form an equation system.
• Solve this system using Singular Value Decomposition (SVD) to find the 3D point in homogeneous coordinates. The resulting point is normalized by dividing by its last coordinate to obtain (X, Y, Z).

4. Color Extraction

For visualization, extract the color information for each 2D point from the first image. Each 2D point corresponds to a pixel in the image, and the RGB color values are stored for each 3D point.

5. Plotting the 3D Points

Finally, plot the 3D points using a scatter plot, where the coordinates (X, Y, Z) are colored according to the extracted pixel color information. This produces a colored 3D point cloud representing the scene.

2. Results

Name	Image
Point Cloud

Q4: Reconstruct your own scene! (20 points)

1. DATASET SOURCE: https://www.agisoft.com/downloads/sample-data/

2D Correspondences	reconstruction