openCV background¶
scratch work for exploring openCV sfm and dense matching pipelines
In [1]:
import numpy as np
import cv2
import os
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#%matplotlib notebook
In [2]:
chess5_db = "./data/chess_5pass/ds1/"
In [3]:
im36 = cv2.imread(os.path.join(chess5_db, "chess_5pass_0036.tif"), 0)
plt.imshow(im36, cmap='gray')
Out[3]:
In [4]:
im37 = cv2.imread(os.path.join(chess5_db, "chess_5pass_0037.tif"), 0)
plt.imshow(im37, cmap='gray')
Out[4]:
In [5]:
im38 = cv2.imread(os.path.join(chess5_db, "chess_5pass_0038.tif"), 0)
plt.imshow(im38, cmap='gray')
Out[5]:
Stereo rectification of pair of images¶
- extract features to match
SIFT features are typically used - but viewing the feature extraction in Alice vision the akaze-liop work better
Alice vision SIFT definition: https://github.com/alicevision/AliceVision/tree/develop/src/aliceVision/feature/sift Alice vision AKAZE definition: https://github.com/alicevision/AliceVision/blob/develop/src/aliceVision/feature/akaze/descriptorLIOP.cpp
A(ccelerated) KAZE feature extraction is an implementation of this with opencv details here
In [6]:
sift = cv2.xfeatures2d.SIFT_create()
akaze = cv2.AKAZE_create()
In [7]:
kp36_sift, des36_sift = sift.detectAndCompute(im36,None)
im36_kp_sift = np.zeros(im36.shape)
im36_kp_sift = cv2.drawKeypoints(im36, kp36_sift, im36_kp_sift, color=[255, 0, 0])
plt.imshow(im36_kp_sift)
Out[7]:
In [8]:
kp37_sift, des37_sift = sift.detectAndCompute(im37,None)
im37_kp_sift = np.zeros(im37.shape)
im37_kp_sift = cv2.drawKeypoints(im37, kp37_sift, im37_kp_sift, color=[255, 0, 0])
plt.imshow(im37_kp_sift)
Out[8]:
In [9]:
kp36_akaze = akaze.detect(im36,None)
im36_kp_akaze = np.zeros(im36.shape)
im36_kp_akaze = cv2.drawKeypoints(im36, kp36_akaze, im36_kp_akaze, color=[0, 255, 0])
plt.imshow(im36_kp_akaze)
Out[9]:
kp - keypoints with some attributes that can be used for thresholding: pt, response, angle, etc.
des - feature descriptors, for SIFT it's a 128 element vector describing gradients weighted by gaussians in a 16x16 neighborhood
In [ ]:
responses = [kp36_sift[i].response for i in range(len(kp36_sift))]
print("response max: ", max(responses))
print("response min: ", min(responses))
plt.hist(responses)
In [10]:
# match sift features between im36 and im37
# using fast library for approximate nearest neightboors
FLANN_INDEX_KDTREE = 0
index_params = dict(algorithm = FLANN_INDEX_KDTREE,
trees = 5
)
search_params = dict(checks = 50)
flann = cv2.FlannBasedMatcher(index_params, search_params)
matches = flann.knnMatch(des36_sift, des37_sift,k=2)
In [11]:
# all matches
matches1to2 = [m[0] for m in matches]
draw_params = dict(matchColor=(0, 0, 255),
singlePointColor = None,
flags=0
)
img_all_matches = cv2.drawMatches(im36, kp36_sift, im37, kp37_sift, matches1to2, None, **draw_params)
plt.figure(figsize = (20,20))
plt.imshow(img_all_matches)
print(len(matches1to2), "matches")
In [12]:
# filter for only good matches
match_thresh = 0.7
# lowe's ratio test
good_match = []
for i, (m, n) in enumerate(matches):
if m.distance < match_thresh*n.distance:
good_match.append(m)
draw_params = dict(matchColor=(0, 0, 255),
singlePointColor = None,
flags = 0
)
img_good_matches = cv2.drawMatches(im36, kp36_sift, im37, kp37_sift, good_match, None, **draw_params)
plt.figure(figsize = (15,15))
plt.imshow(img_good_matches)
print(len(good_match), "matches")
In [13]:
# find essential matrix
im36_matches = np.float32([kp36_sift[m.queryIdx].pt for m in good_match])
im37_matches = np.float32([kp37_sift[m.trainIdx].pt for m in good_match])
F, mask = cv2.findFundamentalMat(im36_matches, im37_matches, cv2.RANSAC, 3, 0.8)
In [14]:
rect, h1, h2 = cv2.stereoRectifyUncalibrated(im36_matches, im37_matches, F, im36.shape[:2])
In [15]:
h1
Out[15]:
In [16]:
h2
Out[16]:
In [17]:
np.array(kp36_sift[good_match[0].queryIdx].pt).astype(int)
Out[17]:
In [18]:
x = range(0, im36.shape[0])
y = range(0, im36.shape[1])
vx, vy = np.meshgrid(x, y)
pt_map = np.stack((vx, vy, np.ones(vx.shape)), axis=2)
In [20]:
homo1 = pt_map@h1
norm1 = np.ones(homo1.shape)
norm1[:, :, 0] = homo1[:, :, 2]
norm1[:, :, 1] = homo1[:, :, 2]
norm1[:, :, 2] = homo1[:, :, 2]
homo_norm1 = np.divide(homo1, norm1)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(pt_map[:,:, 0], pt_map[:, :, 1], pt_map[:, :, 2], color='r')
ax.plot_surface(homo1[:, :, 0], homo1[:, :, 1], homo1[:, :, 2], color='b')
#ax.plot_surface(homo_norm1[:, :, 0], homo_norm1[:, :, 1], homo_norm1[:, :, 2], color='c')
Out[20]:
In [ ]:
homo_norm1
In [21]:
im36_out = cv2.warpPerspective(im36, h1, im36.shape[:2])
im37_out = cv2.warpPerspective(im37, h2, im37.shape[:2])
fig, axes = plt.subplots(1,2, figsize=(15, 10))
axes[0].imshow(im36_out, cmap="gray")
axes[1].imshow(im37_out, cmap='gray')
axes[0].axhline(800)
axes[1].axhline(800)
axes[0].axhline(600)
axes[1].axhline(600)
Out[21]:
In [ ]:
In [27]:
# SGMB
# compute disparity
stereoSGBM = cv2.StereoSGBM_create(minDisparity=0,
numDisparities = 64,
blockSize = 3,
P1 = 8*1*3*3,
P2 = 32*1*3*3
)
dispSGBM = stereoSGBM.compute(im36_out, im37_out).astype(np.float32) / 16
In [28]:
plt.imshow(dispSGBM, 'gray')
plt.colorbar()
Out[28]:
In [32]:
# BM
# computer disparity
stereoBM = cv2.StereoBM_create(numDisparities=16, blockSize=15)
dispBM = stereoBM.compute(im36_out, im37_out)
In [33]:
plt.imshow(dispBM, 'gray')
plt.colorbar()
Out[33]:
In [34]:
dispBM.min()
Out[34]:
In [35]:
# 3D reconstruction
h, w = im36_out.shape
f = 0.8*w # focal length guess
Q = np.float32([[1, 0, 0, -0.5*w],
[0, -1, 0, 0.5*h],
[0, 0, 0, -f],
[0, 0, 1, 0]
])
points = cv2.reprojectImageTo3D(dispSGBM, Q)
In [36]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:, :, 0], points[:, :, 1], points[:, :, 2], color='c')
Out[36]:
In [ ]: