In [18]:
import numpy as np
import cv2
import random
from tqdm import tqdm
import matplotlib.pyplot as plt
import os, pycolmap
import open3d as o3d
from os import mkdir
from pathlib import Path
from IPython.display import display, HTML, Image
In [26]:
####### QUESTION 1
### Functions
def find_F_eightpt(corr1, corr2):
N = corr1.shape[0]
h_corr1 = np.hstack((corr1,np.ones((N,1))))
h_corr2 = np.hstack((corr2,np.ones((N,1))))
x0 = np.mean(h_corr1[:,:2],axis=0)
d_avg = np.sum(np.linalg.norm(h_corr1[:,:2],axis=1))/N
s = np.sqrt(2)/d_avg
T = np.array([[s,0,-s*x0[0]],[0,s,-s*x0[1]],[0,0,1]])
x0 = np.mean(h_corr2[:,:2],axis=0)
d_avg = np.sum(np.linalg.norm(h_corr2[:,:2],axis=1))/N
s = np.sqrt(2)/d_avg
Tpr = np.array([[s,0,-s*x0[0]],[0,s,-s*x0[1]],[0,0,1]])
h_corr1 = (T@h_corr1.T).T
h_corr2 = (Tpr@h_corr2.T).T
A = np.zeros((N,9))
for i in range(N):
x = h_corr1[i,:].reshape((3,1))
xp = h_corr2[i,:]
A[i,:] = np.hstack((xp[0]*x.T, xp[1]*x.T, xp[2]*x.T))
u, s, vt = np.linalg.svd(A)
f = vt[-1,:]/vt[-1,-1]
F = f.reshape((3,3))
u, s, vt = np.linalg.svd(F)
F = u@np.diag([s[0],s[1],0])@vt
return Tpr.T@F@T
def find_E_eightpt(corr1,corr2,K1,K2):
N = corr1.shape[0]
h_corr1 = (np.linalg.inv(K1)@np.hstack((corr1,np.ones((N,1)))).T).T
h_corr2 = (np.linalg.inv(K2)@np.hstack((corr1,np.ones((N,1)))).T).T
A = np.zeros((N,9))
for i in range(N):
x = h_corr1[i,:].reshape((3,1))
xp = h_corr2[i,:]
A[i,:] = np.hstack((xp[0]*x.T, xp[1]*x.T, xp[2]*x.T))
u, s, vt = np.linalg.svd(A)
e = vt[-1,:]/vt[-1,-1]
E = e.reshape((3,3))
u, s, vt = np.linalg.svd(E)
avg_ = (s[0]+s[1])/2
E = u@np.diag([avg_,avg_,0])@vt
return(E)
def find_F_sevenpt(corr1,corr2):
N = corr1.shape[0]
h_corr1 = np.hstack((corr1,np.ones((N,1))))
h_corr2 = np.hstack((corr2,np.ones((N,1))))
A = np.zeros((N,9))
for i in range(N):
x = h_corr1[i,:].reshape((3,1))
xp = h_corr2[i,:]
A[i,:] = np.hstack((xp[0]*x.T, xp[1]*x.T, xp[2]*x.T))
u, s, vt = np.linalg.svd(A)
f1 = vt[-2,:]/vt[-2,-1]
F1 = f1.reshape((3,3))
u, s, vt = np.linalg.svd(F1)
F1 = u@np.diag([s[0],s[1],0])@vt
f1 = F1.reshape((9,1))
u, s, vt = np.linalg.svd(A)
f2 = vt[-1,:]/vt[-1,-1]
F2 = f2.reshape((3,3))
u, s, vt = np.linalg.svd(F2)
F2 = u@np.diag([s[0],s[1],0])@vt
f2 = f2.reshape((9,1))
fc = f1+f2
coeffs = np.ones((4,))
coeffs[0] = fc[0]*(fc[4]*fc[8]-fc[5]*fc[7])-fc[1]*(fc[3]*fc[8]-fc[5]*fc[6])+fc[2]*(fc[3]*fc[7]-fc[4]*fc[6])
coeffs[1] = -np.sum([np.array([f2[0],f2[4],f2[8]]).T@np.array([fc[8]+fc[4],fc[0]+fc[8],fc[0]+fc[4]]),
np.array([f2[0],f2[5],f2[7]]).T@np.array([fc[7]+fc[5],fc[0]+fc[7],fc[0]+fc[5]]),
-np.array([f2[1],f2[3],f2[8]]).T@np.array([fc[3]+fc[8],fc[1]+fc[8],fc[3]+fc[1]]),
-np.array([f2[1],f2[5],f2[6]]).T@np.array([fc[5]+fc[6],fc[1]+fc[5],fc[6]+fc[1]]),
np.array([f2[2],f2[3],f2[7]]).T@np.array([fc[3]+fc[7],fc[2]+fc[7],fc[3]+fc[2]]),
np.array([f2[2],f2[4],f2[6]]).T@np.array([fc[4]+fc[6],fc[2]+fc[4],fc[6]+fc[2]])])
coeffs[2] = np.sum([np.array([f2[4]*f2[8],f2[0]*f2[8],f2[0]*f2[4]]).T@np.array([fc[0],fc[4],fc[8]]),
np.array([f2[5]*f2[7],f2[0]*f2[7],f2[0]*f2[5]]).T@np.array([fc[0],fc[5],fc[7]]),
-np.array([f2[3]*f2[8],f2[1]*f2[8],f2[1]*f2[3]]).T@np.array([fc[1],fc[3],fc[8]]),
-np.array([f2[5]*f2[6],f2[1]*f2[6],f2[1]*f2[5]]).T@np.array([fc[1],fc[5],fc[6]]),
np.array([f2[3]*f2[7],f2[2]*f2[7],f2[2]*f2[3]]).T@np.array([fc[2],fc[3],fc[7]]),
np.array([f2[4]*f2[6],f2[2]*f2[6],f2[2]*f2[4]]).T@np.array([fc[2],fc[4],fc[6]])])
coeffs[3] = -f2[0]*(f2[4]*f2[8]-f2[5]*f2[7])+f2[1]*(f2[3]*f2[8]-f2[5]*f2[6])-f2[2]*(f2[4]*f2[7]-f2[4]*f2[6])
soln = np.roots(coeffs)
F_s = [soln[i]*F1+(1-soln[i])*F2 for i in range(len(soln))]
F_s = [F/F[-1,-1] for F in F_s]
F_s = []
for sol in soln:
F = sol*F1+(1-sol)*F2
F = F/F[-1,-1]
if (h_corr2[0].T@F@h_corr1[0])**2 < 1 and not np.iscomplexobj(F):
F_s.append(F)
try:
errors = np.zeros((len(F_s),1))
for j in range(N):
for i, F in enumerate(F_s):
errors[i] += (h_corr2[j].T@F@h_corr1[j])**2
F = F_s[np.argmin(errors)]
u, s, vt = np.linalg.svd(F)
F = u@np.diag([s[0],s[1],0])@vt
return F
except:
return np.ones((3,3))
def return_plots(im1,im2,corr1,corr2,F):
def draw_epipolar_lines(image, F, lines):
N = lines.shape[0]
colors = [(255,0,0),(0,255,0),(0,0,255),(255,255,0),(255,0,255),(0,255,255),(0,0,0),(255,255,255)]
for i in range(min(N,len(colors))):
l = lines[i]
# Define the range of x' values (for the image width)
x_prime = np.array([0, image.shape[1] - 1])
# Calculate the corresponding y' values
y_prime = -(l[0] * x_prime + l[2]) / l[1]
# Convert to integer coordinates for drawing
pt1 = (int(x_prime[0]), int(y_prime[0]))
pt2 = (int(x_prime[1]), int(y_prime[1]))
# Draw the epipolar line
cv2.line(image, pt1, pt2, color=colors[i], thickness=7) # Red line
return image
def draw_epipolar_points(image, points):
N = points.shape[0]
colors = [(255,0,0),(0,255,0),(0,0,255),(255,255,0),(255,0,255),(0,255,255),(0,0,0),(255,255,255)]
for i in range(min(N,len(colors))):
cv2.circle(image, tuple(points[i].astype(np.int64)), 10, colors[i], -1)
return image
N = corr1.shape[0]
h_corr1 = np.hstack((corr1,np.ones((N,1))))
h_corr2 = np.hstack((corr2,np.ones((N,1))))
lines = np.ones((N,3))
for i in range(N):
line = F@h_corr1[i].T
line /= line[-1]
lines[i] = line
im1 = draw_epipolar_points(im1, corr1)
im2 = draw_epipolar_lines(im2, F, lines)
plt.figure(figsize=(10, 10), dpi=80)
plt.subplot(1,2,1)
plt.imshow(cv2.cvtColor(im1, cv2.COLOR_BGR2RGB))
plt.title("Image 1 - Points")
plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(cv2.cvtColor(im2, cv2.COLOR_BGR2RGB))
plt.title("Image 2 - Lines")
plt.show()
In [27]:
corr_file = "./data/q1a/bench/corresp.npz"
intrinsics_file = "./data/q1a/bench/intrinsics.npz"
im1_file = "./data/q1a/bench/image1.jpg"
im2_file = "./data/q1a/bench/image2.jpg"
corresp = np.load(corr_file)
intrin = np.load(intrinsics_file)
im1 = cv2.imread(im1_file)
im2 = cv2.imread(im2_file)
corr1 = corresp['pts1']
corr2 = corresp['pts2']
K1 = intrin['K1']
K2 = intrin['K2']
Explanation - Problem Q1A1¶
In order to find the eight-point reconstruction of the fundamental matrix, I made the image coordinates homogeneous and normalized, as well as creating an A matrix that followed the linear constraint (x')^T * F * x = 0. With the A matrix constructed with eight necessary points in order to find a solution, the A matrix was passed through an SVD and a flattened f vector was found from the last singular vector. Lastly, rank-F was found using a rank-2 approximation (SVD(F) -> U * diag(d_1, d_2, 0) * V^T = F. F was "un-normalized" and was passed for plotting.
In [28]:
### A1
F = find_F_eightpt(corr1, corr2)
print(" --------- Q1A (Fundamental) ---------")
print(F)
return_plots(im1, im2, corr1, corr2, F)
--------- Q1A (Fundamental) --------- [[ 1.50229865e-05 3.19006224e-04 -2.61477507e-02] [ 4.99396287e-05 -2.86055451e-05 -3.13479290e+00] [ 1.34641936e-01 2.89778705e+00 -1.25962199e+02]]
Explanation - Problem Q1A2¶
In order to find the eight-point reconstruction of the essential matrix, the same process was followed as the above fundamental matrix, with the incorporation of the relevant intrinsics matrices. Each of the points was made homogeneous and multiplied by the intrinsics matrix for normalization. Then, the rest of the F eight-point algorithm was followed, with the only difference being the rank 2 approximation, with E = U * diag(d_avg, d_avg, 0) * V^T. E was passed for plotting.
In [29]:
### A2
E = find_E_eightpt(corr1, corr2, K1, K2)
print(" --------- Q1A (Essential) ---------")
print(E)
--------- Q1A (Essential) --------- [[-5.77317706e-01 -2.90883845e+03 2.49343487e+03] [ 2.90971269e+03 8.62381999e-01 2.23890104e+03] [-2.49220018e+03 -2.23829395e+03 7.14935646e-01]]
Explanation - Problem Q1B¶
For the seven-point algorithm, the A matrix was constructed similarly (with x'.T * F * x = 0 as the basis constraint), and SVD was used to decompose A. From the decomposition, the two singular vectors of A were taken for finding the root that solved det(lambdaF1 + (1-lambda)F2) = 0. The particular polynomial coefficient were manually calculated and made into a 4-element vector, which was passed through np.roots for possible solutions to the polynomial. Imaginary roots were excluded, and the F that had the least error for the points used for correspondence was returned as the fundamental matrix for the set of seven points.
In [30]:
### B
corr_file = "./data/q1b/hydrant/corresp.npz"
im1_file = "./data/q1b/hydrant/image1.jpg"
im2_file = "./data/q1b/hydrant/image2.jpg"
corresp = np.load(corr_file)
im1 = cv2.imread(im1_file)
im2 = cv2.imread(im2_file)
corr1 = corresp['pts1']
corr2 = corresp['pts2']
F = find_F_sevenpt(corr1, corr2)
return_plots(im1, im2, corr1, corr2, F)
print(" --------- Q1B ---------")
print(F)
--------- Q1B --------- [[ 9.12671347e-07 7.72710923e-06 -4.66445103e-03] [-4.72074193e-06 1.11051900e-06 -5.62393695e-03] [ 3.48580503e-03 3.53316889e-03 1.00000000e+00]]
In [31]:
####### QUESTION 2
### Functions
def ransacF(corr_file, n, t, iter_):
def inlier_find(corr1, corr2, F, t):
N = corr1.shape[0]
count = 0
for i in range(N):
line = F @ corr1[i]
err = np.abs(corr2[i].T @ line)
if err < t:
count +=1
if count/N > .4:
return 1,count
else:
return 0,count
n = 7
t = 1e-2
corresp = np.load(corr_file)
im1 = cv2.imread(im1_file)
im2 = cv2.imread(im2_file)
corr1 = corresp['pts1']
corr2 = corresp['pts2']
N = corr1.shape[0]
h_corr1 = np.hstack((corr1,np.ones((N,1))))
h_corr2 = np.hstack((corr2,np.ones((N,1))))
list_N = list(range(N))
best_inliers = 0
best_F = np.eye(3)
inl_list = np.zeros((iter_,1))
print(" --------- RANSAC ---------")
for i in tqdm(range(iter_)):
if n == 8:
random.shuffle(list_N)
idxs = list_N[:8]
corr1_ = corr1[idxs]
corr2_ = corr2[idxs]
F = find_F_eightpt(corr1_, corr2_)
inliers = inlier_find(h_corr1,h_corr2,F, t)[1]
if inliers > best_inliers:
best_F = F[:,:]
best_inliers = inliers
elif n == 7:
random.shuffle(list_N)
idxs = list_N[:7]
corr1_ = corr1[idxs]
corr2_ = corr2[idxs]
F = find_F_sevenpt(corr1_, corr2_)
inliers = inlier_find(h_corr1,h_corr2,F, t)[1]
if not np.all(F == np.ones((3,3))):
if inliers > best_inliers:
best_F = F[:,:]
best_inliers = inliers
inl_list[i] = best_inliers/N
F_ransac = best_F.reshape((3,3))
u, s, vt = np.linalg.svd(F_ransac)
F_ransac = u@np.diag([s[0],s[1],0])@vt
return F_ransac, inl_list
def return_plots_2(corr_file, t, iter_7, iter_8):
F8, inliers_8 = ransacF(corr_file, 8, t, iter_8)
F7, inliers_7 = ransacF(corr_file, 7, t, iter_7)
print(F)
plt.figure()
plt.plot(np.array(list(range(iter_7))),inliers_7,label="7-point")
plt.plot(np.array(list(range(iter_8))),inliers_8,label='8-point')
plt.legend()
plt.title("Inlier Prediction Updates")
plt.ylabel("Percent Inliers")
plt.xlabel("Iterations")
plt.show()
return_plots(im1, im2, corr1, corr2, F8)
return None
Explanation - Problem Q2¶
For the RANSAC algorithm, the minimum number of points (7 or 8 depending on the algorithm) were randomly sampled from the test set and used for inlier finding. Once a set of points came up with a fundamental matrix, it was tested against the entire set of points to see if the distance to the epipolar line |x_i.T @ l_i| was less than some threshold (0.01). The count of inliers was maintained, and the fundamental matrix that was returned at the end of the function was the F with the most inliers. The max number of iteration was a probabilistically determined number that found the number of iterations necessary for a certain probability of all inliers being in one set (probability used was 0.01). The F_ransac was returned and used for plotting.
In [42]:
## 2
corr_file = "./data/q1b/hydrant/corresp_noisy.npz"
im1_file = "./data/q1b/hydrant/image1.jpg"
im2_file = "./data/q1b/hydrant/image2.jpg"
corr_smooth_file = "./data/q1b/hydrant/corresp.npz"
corresp = np.load(corr_smooth_file)
im1 = cv2.imread(im1_file)
im2 = cv2.imread(im2_file)
corr1 = corresp['pts1']
corr2 = corresp['pts2']
n = 7
t = 0.1
iter_7 = int(np.log(.05)/np.log((1-(.4)**n)))
iter_8 = int(np.log(.05)/np.log((1-(.4)**8)))
F, _ = ransacF(corr_file, 7, t, iter_7)
print(" --------- Q2A ---------")
print(F)
return_plots_2(corr_file, t, iter_7, iter_8)
--------- RANSAC ---------
100%|███████████████████████████████████████| 1826/1826 [00:25<00:00, 72.17it/s]
--------- Q2A --------- [[ 6.14364703e-07 4.14631771e-06 -5.61940705e-03] [-3.44283625e-06 1.02320809e-06 -2.01847395e-03] [ 5.26150334e-03 -5.44851249e-05 1.00000000e+00]] --------- RANSAC ---------
100%|███████████████████████████████████████| 4569/4569 [01:02<00:00, 73.28it/s]
--------- RANSAC ---------
100%|███████████████████████████████████████| 1826/1826 [00:25<00:00, 72.79it/s]
[[ 6.14364703e-07 4.14631771e-06 -5.61940705e-03] [-3.44283625e-06 1.02320809e-06 -2.01847395e-03] [ 5.26150334e-03 -5.44851249e-05 1.00000000e+00]]
In [33]:
####### QUESTION 3 - TRIANGULATION
### Functions
def triangulation(P1_file,P2_files, corr1_file, corr2_file):
def cross_matrix(x):
return np.array([[0,-x[2],x[1]],[x[2],0,-x[0]],[-x[1],x[0],0]])
P1 = np.load(P1_file)
P2 = np.load(P2_file)
corr1 = np.load(corr1_file)
corr2 = np.load(corr2_file)
N = corr1.shape[0]
h_corr1 = np.hstack((corr1,np.ones((N,1))))
h_corr2 = np.hstack((corr2,np.ones((N,1))))
print(" --------- TRIANGULATION ---------")
X = np.ones((N,4))
for i in tqdm(range(N)):
A = np.vstack(((cross_matrix(h_corr1[i])@P1)[:2,:],(cross_matrix(h_corr2[i])@P2)[:2,:]))
_,_,vt = np.linalg.svd(A)
X[i,:] = vt[-1,:]/vt[-1,-1]
return X
def return_plots3(corr1_file, corr2_file, im1_file, im2_file, X):
corr1 = np.load(corr1_file)
im1 = cv2.imread(im1_file)
corr2 = np.load(corr2_file)
im2 = cv2.imread(im2_file)
N = corr1.shape[0]
colors = np.zeros((N,3))
for i in range(N):
try:
color = (im1[corr1[i][1],corr1[i][0]])/255
colors[i] = [color[2],color[1],color[0]]
except:
color = (im2[corr2[i][1],corr2[i][0]])/255
colors[i] = [color[2],color[1],color[0]]
# Create a figure for plotting
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot the sampled points
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=colors, marker='o')
# Label axes
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
ax.view_init(155, 90)
# Show the plot
plt.show(block=True)
Explanation - Problem Q3¶
For triangulation, the basis for finding the 3d points was the formation of a matrix A for singular value decomposition. Following the constraint that X cross Px = 0, the A was constructed from the first two rows of [x]_x * P and [x']_x * P' vertically stacked. The singular vector of A was the 3D coordinate X, which was normalized based on its homogenous coordinate. This was repeated for all points X and the colors were maintained for the final plotting.
In [34]:
im1_file = "./data/q3/img1.jpg"
im2_file = "./data/q3/img2.jpg"
P1_file = "./data/q3/P1.npy"
P2_file = "./data/q3/P2.npy"
corr1_file = "./data/q3/pts1.npy"
corr2_file = "./data/q3/pts2.npy"
X = triangulation(P1_file,P2_file,corr1_file, corr2_file)
return_plots3(corr1_file, corr2_file, im1_file, im2_file, X)
--------- TRIANGULATION ---------
100%|████████████████████████████████| 100000/100000 [00:06<00:00, 16633.90it/s]
Explanation - Problem Q4¶
For the sparse reconstruction from a custom scene, I followed the tutorial given by the PyCOLMAP Github page, which describes setting up the requisite data paths with pathlib. From there, you extract features, match, and then map incrementally to find the basic reconstruction. The map was then output and reformatting into a PLY file, which were imported into MeshLab for viewing and .gif generation.
In [35]:
####### QUESTION 4 - RECONSTRUCT SCENE
test = "full_mug"
output_path = Path.cwd() / ("data/q4_scenes/"+test+"/tests")
output_path.mkdir(exist_ok=True)
im_dir = Path.cwd() / ("data/q4_scenes/"+test+"/pics")
database_path = output_path / (test+".db")
pycolmap.extract_features(database_path, im_dir)
pycolmap.match_exhaustive(database_path)
maps = pycolmap.incremental_mapping(database_path, im_dir, output_path)
maps[0].write(output_path)
reconstruction = pycolmap.Reconstruction(str(output_path))
print(reconstruction.summary())
# reconstruction.export_PLY("data/q4_scenes/"+test+"/"+test+".ply") # PLY format
W20241031 21:35:39.270212 0x7ff84d564b80 feature_extraction.cc:406] Your current options use the maximum number of threads on the machine to extract features. Extracting SIFT features on the CPU can consume a lot of RAM per thread for large images. Consider reducing the maximum image size and/or the first octave or manually limit the number of extraction threads. Ignore this warning, if your machine has sufficient memory for the current settings. I20241031 21:35:39.272342 0x70000d783000 misc.cc:198] ============================================================================== Feature extraction ============================================================================== I20241031 21:35:39.272621 0x7000106b9000 sift.cc:722] Creating SIFT CPU feature extractor I20241031 21:35:39.272648 0x700012848000 sift.cc:722] Creating SIFT CPU feature extractor I20241031 21:35:39.272751 0x7000128cb000 sift.cc:722] Creating SIFT CPU feature extractor I20241031 21:35:39.272843 0x70001294e000 sift.cc:722] Creating SIFT CPU feature extractor I20241031 21:35:39.272863 0x700012a54000 sift.cc:722] Creating SIFT CPU feature extractor I20241031 21:35:39.272869 0x700012ad7000 sift.cc:722] Creating SIFT CPU feature extractor I20241031 21:35:39.272879 0x700012b5a000 sift.cc:722] Creating SIFT CPU feature extractor I20241031 21:35:39.272912 0x7000129d1000 sift.cc:722] Creating SIFT CPU feature extractor I20241031 21:35:39.276499 0x700012bdd000 feature_extraction.cc:257] Processed file [1/13] I20241031 21:35:39.276556 0x700012bdd000 feature_extraction.cc:260] Name: .DS_Store E20241031 21:35:39.276566 0x700012bdd000 feature_extraction.cc:266] Failed to read image file format. I20241031 21:35:39.279311 0x700012bdd000 feature_extraction.cc:257] Processed file [2/13] I20241031 21:35:39.279371 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1009.png I20241031 21:35:39.279385 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.279394 0x700012bdd000 feature_extraction.cc:257] Processed file [3/13] I20241031 21:35:39.279402 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1010.png I20241031 21:35:39.279409 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.280218 0x700012bdd000 feature_extraction.cc:257] Processed file [4/13] I20241031 21:35:39.280349 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1011.png I20241031 21:35:39.280374 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.280577 0x700012bdd000 feature_extraction.cc:257] Processed file [5/13] I20241031 21:35:39.280596 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1012.png I20241031 21:35:39.280607 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.281121 0x700012bdd000 feature_extraction.cc:257] Processed file [6/13] I20241031 21:35:39.281152 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1013.png I20241031 21:35:39.281165 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.281826 0x700012bdd000 feature_extraction.cc:257] Processed file [7/13] I20241031 21:35:39.281855 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1014.png I20241031 21:35:39.281868 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.281877 0x700012bdd000 feature_extraction.cc:257] Processed file [8/13] I20241031 21:35:39.281885 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1015.png I20241031 21:35:39.281894 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.282231 0x700012bdd000 feature_extraction.cc:257] Processed file [9/13] I20241031 21:35:39.282254 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1016.png I20241031 21:35:39.282265 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.283453 0x700012bdd000 feature_extraction.cc:257] Processed file [10/13] I20241031 21:35:39.283502 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1017.png I20241031 21:35:39.283521 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.284168 0x700012bdd000 feature_extraction.cc:257] Processed file [11/13] I20241031 21:35:39.284206 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1018.png I20241031 21:35:39.284320 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.284332 0x700012bdd000 feature_extraction.cc:257] Processed file [12/13] I20241031 21:35:39.284341 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1019.png I20241031 21:35:39.284350 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.285024 0x700012bdd000 feature_extraction.cc:257] Processed file [13/13] I20241031 21:35:39.285053 0x700012bdd000 feature_extraction.cc:260] Name: IMG_1020.png I20241031 21:35:39.285065 0x700012bdd000 feature_extraction.cc:264] SKIP: Features for image already extracted. I20241031 21:35:39.294452 0x70000d783000 timer.cc:91] Elapsed time: 0.000 [minutes] I20241031 21:35:39.317381 0x70000d783000 misc.cc:198] ============================================================================== Feature matching ============================================================================== I20241031 21:35:39.317671 0x7000102a1000 sift.cc:1457] Creating SIFT CPU feature matcher I20241031 21:35:39.317711 0x7000103a7000 sift.cc:1457] Creating SIFT CPU feature matcher I20241031 21:35:39.317719 0x70001042a000 sift.cc:1457] Creating SIFT CPU feature matcher I20241031 21:35:39.317744 0x7000104ad000 sift.cc:1457] Creating SIFT CPU feature matcher I20241031 21:35:39.317711 0x700010324000 sift.cc:1457] Creating SIFT CPU feature matcher I20241031 21:35:39.317782 0x700010530000 sift.cc:1457] Creating SIFT CPU feature matcher I20241031 21:35:39.317826 0x7000105b3000 sift.cc:1457] Creating SIFT CPU feature matcher I20241031 21:35:39.317863 0x700010636000 sift.cc:1457] Creating SIFT CPU feature matcher I20241031 21:35:39.319404 0x70000d783000 pairing.cc:168] Generating exhaustive image pairs... I20241031 21:35:39.319464 0x70000d783000 pairing.cc:202] Matching block [1/1, 1/1] I20241031 21:35:39.326544 0x70000d783000 feature_matching.cc:46] in 0.007s I20241031 21:35:39.326767 0x70000d783000 timer.cc:91] Elapsed time: 0.000 [minutes] I20241031 21:35:39.331086 0x7ff84d564b80 incremental_mapper.cc:225] Loading database I20241031 21:35:39.332702 0x7ff84d564b80 database_cache.cc:65] Loading cameras... I20241031 21:35:39.332799 0x7ff84d564b80 database_cache.cc:75] 22 in 0.000s I20241031 21:35:39.332842 0x7ff84d564b80 database_cache.cc:83] Loading matches... I20241031 21:35:39.333053 0x7ff84d564b80 database_cache.cc:89] 45 in 0.000s I20241031 21:35:39.333065 0x7ff84d564b80 database_cache.cc:105] Loading images... I20241031 21:35:39.339338 0x7ff84d564b80 database_cache.cc:155] 12 in 0.006s (connected 12) I20241031 21:35:39.339368 0x7ff84d564b80 database_cache.cc:166] Building correspondence graph... I20241031 21:35:39.342060 0x7ff84d564b80 database_cache.cc:195] in 0.003s (ignored 0) I20241031 21:35:39.342092 0x7ff84d564b80 timer.cc:91] Elapsed time: 0.000 [minutes] I20241031 21:35:39.344340 0x7ff84d564b80 incremental_mapper.cc:265] Finding good initial image pair I20241031 21:35:39.363484 0x7ff84d564b80 incremental_mapper.cc:289] Initializing with image pair #4 and #2 I20241031 21:35:39.364279 0x7ff84d564b80 incremental_mapper.cc:294] Global bundle adjustment I20241031 21:35:39.491224 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #12 (3) I20241031 21:35:39.491261 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 117 / 466 points I20241031 21:35:39.625636 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:39.743803 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #3 (4) I20241031 21:35:39.743908 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 180 / 349 points I20241031 21:35:39.832332 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:39.891758 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #10 (5) I20241031 21:35:39.891794 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 223 / 467 points I20241031 21:35:39.946978 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:40.024879 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #1 (6) I20241031 21:35:40.024919 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 242 / 432 points I20241031 21:35:40.104000 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:40.199806 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #9 (7) I20241031 21:35:40.199851 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 158 / 221 points I20241031 21:35:40.252955 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:40.350378 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #5 (8) I20241031 21:35:40.350426 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 161 / 427 points I20241031 21:35:40.424912 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:40.546217 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #11 (9) I20241031 21:35:40.546257 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 167 / 393 points I20241031 21:35:40.664726 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:40.790471 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #7 (10) I20241031 21:35:40.790511 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 97 / 157 points I20241031 21:35:40.840397 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:40.953651 0x7ff84d564b80 incremental_mapper.cc:375] Registering image #8 (11) I20241031 21:35:40.953695 0x7ff84d564b80 incremental_mapper.cc:378] => Image sees 86 / 131 points I20241031 21:35:40.986616 0x7ff84d564b80 incremental_mapper.cc:41] Retriangulation and Global bundle adjustment I20241031 21:35:41.143079 0x7ff84d564b80 timer.cc:91] Elapsed time: 0.030 [minutes]
Reconstruction: num_reg_images = 11 num_cameras = 11 num_points3D = 632 num_observations = 2800 mean_track_length = 4.43038 mean_observations_per_image = 254.545 mean_reprojection_error = 0.484176
In [21]:
gifPath = Path("/Users/morganmayborne/Downloads/proj3/bobblehead.gif")
# Display GIF in Jupyter, CoLab, IPython
with open(gifPath,'rb') as f:
display(Image(data=f.read(), format='png'))
In [22]:
gifPath = Path("/Users/morganmayborne/Downloads/proj3/full_mug_video.gif")
# Display GIF in Jupyter, CoLab, IPython
with open(gifPath,'rb') as f:
display(Image(data=f.read(), format='png'))
Explanation - Problem Q5¶
With the initial code given on the problem set guide, kp and des were computed for both images given. From here, I chose to use a Brute-Force matcher to map the descriptors to image points. I also chose to sort the matched points by distance error, to ensure the top points in the set were stronger matches. Lastly, the points were extracted, and then I ran my Q1 eight_point algorithm on the points to find F and plot the results.
In [38]:
####### QUESTION 5 - AUTOMATIC FUNDAMENTAL MATRIX
# Loading the first image
img1 = cv2.imread("data/q4_scenes/bobblehead/pics/IMG_1021.png")
print(img1.shape)
gray1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)
# Applying SIFT detector to the first image
sift = cv2.xfeatures2d.SIFT_create()
kp1, des1 = sift.detectAndCompute(gray1, None)
# Loading the second image
img2 = cv2.imread("data/q4_scenes/bobblehead/pics/IMG_1022.png")
gray2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)
# Applying SIFT detector to the second image
kp2, des2 = sift.detectAndCompute(gray2, None)
# Creating a BFMatcher object
bf = cv2.BFMatcher(cv2.NORM_L2, crossCheck=True)
# Matching descriptors
matches = bf.match(des1, des2)
# Sorting matches by distance
matches = sorted(matches, key=lambda x: x.distance)
# Extracting the matched keypoints' pixel coordinates
points1 = np.zeros((len(matches), 2), dtype=np.float32)
points2 = np.zeros((len(matches), 2), dtype=np.float32)
for i, match in enumerate(matches):
points1[i, :] = kp1[match.queryIdx].pt
points2[i, :] = kp2[match.trainIdx].pt
print(points1.shape, points2.shape)
F = find_F_eightpt(points1[:8], points2[:8])
print(F)
img1 = cv2.imread("data/q4_scenes/bobblehead/pics/IMG_1021.png")
img2 = cv2.imread("data/q4_scenes/bobblehead/pics/IMG_1022.png")
return_plots(img1, img2, points1, points2, F)
(640, 480, 3) (236, 2) (236, 2) [[-2.06409436e-04 4.86591742e-04 -1.73064853e-01] [-3.13037517e-04 -8.41481194e-05 2.19693489e-03] [ 1.64146111e-01 1.36964923e-02 8.64379957e+00]]
