import numpy as np
import matplotlib.pyplot as plt
import cv2
%matplotlib inline
# Lecture et affichage des deux images :
I12 = plt.imread('0012.ppm')
I15 = plt.imread('0015.ppm')
plt.figure(figsize=(14,7))
plt.subplot(1,2,1)
plt.imshow(I12)
plt.subplot(1,2,2)
plt.imshow(I15)
plt.tight_layout()
plt.show()
# Sélection manuelle de points en correspondance :
# N.B. Ne semble pas fonctionner dans un notebook
plt.figure(figsize=(14,7))
plt.subplot(1,2,1)
plt.imshow(I12)
p12 = plt.ginput() # décommenter pour effectuer la sélection
plt.subplot(1,2,2)
plt.imshow(I15)
p15 = plt.ginput() # décommenter pour effectuer la sélection
plt.tight_layout()
plt.show()
# Deux listes de points en coorespondance :
p12 = np.array([[411, 352, 1],
[472, 423, 1],
[708, 237, 1],
[304, 181, 1],
[520, 654, 1],
[439, 636, 1],
[567, 559, 1],
[667, 516, 1]])
p15 = np.array([[249, 343, 1],
[363, 417, 1],
[382, 261, 1],
[399, 181, 1],
[354, 643, 1],
[387, 615, 1],
[550, 573, 1],
[780, 582, 1]])
# indices de connectivité (pour afficher des segments reliant les points)
c = np.array([[0,1],[1,2],[1,3],[1,4],[1,5],[1,6],[6,7]])
def plotSegments2D(p,c):
v0 = p[c[:,0],:]
v1 = p[c[:,1],:]
plt.plot(np.array([v0[:,0],v1[:,0]]),np.array([v0[:,1],v1[:,1]]),'b')
# Affichage des images et points sélectionnés
plt.figure(figsize=(14,7))
plt.subplot(1,2,1)
plt.imshow(I12)
plt.plot(p12[:,0],p12[:,1],'ob')
plotSegments2D(p12,c)
plt.subplot(1,2,2)
plt.imshow(I15)
plt.plot(p15[:,0],p15[:,1],'o')
plotSegments2D(p15,c)
plt.tight_layout()
plt.show()
F, mask = cv2.findFundamentalMat(p12,p15,cv2.NORM_L2)
F
array([[ 5.15412392e-07, 5.07422627e-06, -1.86351184e-03], [ 3.48057985e-06, 3.16401476e-07, -1.25618615e-02], [-2.78567853e-04, 8.12745149e-03, 1.00000000e+00]])
# Calcul et affichage des droites épipolaires
def drawEpilines(F,points,img):
x = np.array([0,img.shape[1]-1])
plt.imshow(img)
for p in points:
a, b, c = F@p
y = (-a*x-c)/b
plt.plot(x,y)
plt.figure(figsize=(14,7))
drawEpilines(F,p12,I15)
plt.plot(p15[:,0],p15[:,1],'o')
plt.tight_layout()
plt.show()
plt.figure(figsize=(14,7))
drawEpilines(F.T,p15,I12)
plt.plot(p12[:,0],p12[:,1],'o')
plt.tight_layout()
plt.show()
mask
array([[1], [1], [1], [1], [1], [1], [0], [1]], dtype=uint8)
# sélection d'autres points dans la première image
q12 = np.array([[340, 315, 1],
[619, 346, 1],
[460, 343, 1],
[501, 552, 1]
])
plt.figure(figsize=(14,7))
plt.imshow(I12)
plt.plot(q12[:,0],q12[:,1],'ob')
plt.show()
plt.figure(figsize=(14,7))
plt.subplot(1,2,1)
plt.imshow(I12)
plt.plot(q12[:,0],q12[:,1],'ob')
plt.subplot(1,2,2)
plt.imshow(I15)
drawEpilines(F,q12,I15)
plt.tight_layout()
plt.show()
# Matrices de projection
M12 = np.array([[-2999.61,-588.708,-66.03,37412.1],
[-43.6705,185.125,-3044.78,21258.2 ],
[0.0139967,-0.98101,-0.193453,62.6864]])
M15 = np.array([[-1056.3,-2876.12,-172.27,32219.2],
[-120.743,142.125,-3045.28,22434.5 ],
[0.846132,-0.494588,-0.198602,62.0932]])
# Reconstruction 3D des points
def crossmatrix(p):
return np.array([[0,-p[2],p[1]],[p[2],0,-p[0]],[-p[1],p[0],0]])
def plotSegments3D(p,c,ax):
for i in range(len(c)):
ax.plot([p[c[i,0],0],p[c[i,1],0]],[p[c[i,0],1],p[c[i,1],1]],[p[c[i,0],2],p[c[i,1],2]],'b')
Xrec = np.zeros((len(p12),3))
for i in range(len(p12)):
p = p12[i]
q = p15[i]
M = np.vstack((crossmatrix(p)@M12,crossmatrix(q)@M15))
A = M[:,:3]
b = M[:,3][:,None]
X, _, _, _ = np.linalg.lstsq(A,-b, rcond=None)
Xrec[i,:] = X.flatten()
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111, projection='3d')
plotSegments3D(Xrec,c,ax)
ax.plot(Xrec[:,0],Xrec[:,1],Xrec[:,2],'o')
plt.tight_layout()
plt.show()