import numpy as np
import matplotlib.pyplot as plt
import cv2
%matplotlib notebook
###########################################
# Reprise des éléments de la partie 1
###########################################
# Fonctions utiles
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 # droite épipolaire : équation ax+by+c=0
y = (-a*x-c)/b
plt.plot(x,y)
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')
def Reconstr3D(points_a, points_b, Ma, Mb):
# calcul des coordonnées 3D des points connaissant :
# points_a, points_b : matrices (n,3) des coordonnées des points dans l'image a et l'image b
# (N.B. coordonnées homogènes sous la forme (u,v,1))
# Ma, Mb : matrices (3,4) de projection des deux caméras
# Renvoie Xrec matrice (n,3) des coordonnées 3D
Xrec = np.zeros((len(points_a),3))
for i in range(len(points_b)):
p = points_a[i]
q = points_b[i]
M = np.vstack((crossmatrix(p)@Ma,crossmatrix(q)@Mb))
A = M[:,:3]
b = M[:,3][:,None]
X, _, _, _ = np.linalg.lstsq(A,-b, rcond=None)
Xrec[i,:] = X.flatten()
return Xrec
# Lecture des deux images
Ia = plt.imread('0012.ppm')
Ib = plt.imread('0015.ppm')
# Matrices de projection
Ma = 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]])
Mb = 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]])
# points en correspondance repérés manuellement
pa = 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]])
pb = 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]])
# Estimation de la matrice fondamentale
F, _ = cv2.findFundamentalMat(pa,pb,cv2.FM_LMEDS)
# Question 1 : Calculer et afficher les points SIFT des deux images,
# puis calculer et afficher des points en correspondance entre les deux images pour les descripteurs
# Création d'un objet de la classe SIFT
sift = cv2.SIFT_create()
# Détection des points SIFT sur les deux images
kpa = sift.detect(Ia,None)
kpb = sift.detect(Ib,None)
# Affichage des images avec leurs points
resa = Ia.copy()
resb = Ib.copy()
resa = cv2.drawKeypoints(Ia,kpa,resa,flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
resb = cv2.drawKeypoints(Ib,kpb,resb,flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
%matplotlib inline
plt.figure(figsize=(16,8))
plt.subplot(1,2,1)
plt.imshow(resa)
plt.subplot(1,2,2)
plt.imshow(resb)
plt.tight_layout()
plt.show()
# Calcul des descripteurs SIFT pour les deux images
kpa, dsca = sift.compute(Ia,kpa)
kpb, dscb = sift.compute(Ib,kpb)
# Création d'un objet pour la mise en correspondance des keypoints
bf = cv2.BFMatcher(cv2.NORM_L2)
# Calcul des correspondances
matches = bf.match(dsca, dscb)
# Affichage des correspondances
image_matches = cv2.drawMatches(Ia, kpa, Ib, kpb, matches, None)
plt.figure(figsize=(16,8))
plt.imshow(image_matches)
plt.tight_layout()
plt.show()
# Question 2 : Mêmes questions en remplaçant la méthode SIFT par la méthode ORB
# Création d'un objet de la classe ORB
orb = cv2.ORB_create()
# Détection des points ORB sur les deux images
kpa = orb.detect(Ia,None)
kpb = orb.detect(Ib,None)
# Affichage des images avec leurs points
resa = Ia.copy()
resb = Ib.copy()
resa = cv2.drawKeypoints(Ia,kpa,resa,flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
resb = cv2.drawKeypoints(Ib,kpb,resb,flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
plt.figure(figsize=(16,8))
plt.subplot(1,2,1)
plt.imshow(resa)
plt.subplot(1,2,2)
plt.imshow(resb)
plt.tight_layout()
plt.show()
# Calcul des descripteurs ORB pour les deux images
kpa, dsca = orb.compute(Ia,kpa)
kpb, dscb = orb.compute(Ib,kpb)
# Création d'un objet pour la mise en correspondance des keypoints
bf = cv2.BFMatcher(cv2.NORM_L2)
# Calcul des correspondances
matches = bf.match(dsca, dscb)
# Affichage des correspondances
image_matches = cv2.drawMatches(Ia, kpa, Ib, kpb, matches, None)
plt.figure(figsize=(16,8))
plt.imshow(image_matches)
plt.tight_layout()
plt.show()
# conversion de kpa et kpb en NumPy array de taille (1266,3) contenant les coordonnées (u,v,1) de chaque keypoint
pa = np.array([np.array(kpa[i].pt+(1,)) for i in range(len(kpa))])
pb = np.array([np.array(kpb[i].pt+(1,)) for i in range(len(kpb))])
# Question 3 : Afficher sur la 2e image les droites épipolaires correspondant aux points ORB de la 1e image,
# et inversement.
plt.figure(figsize=(14,7))
drawEpilines(F,pa,Ib)
plt.plot(pb[:,0],pb[:,1],'o')
plt.tight_layout()
plt.show()
plt.figure(figsize=(14,7))
drawEpilines(F.T,pb,Ia)
plt.plot(pa[:,0],pa[:,1],'o')
plt.tight_layout()
plt.show()
# Question 7 : Calculer une estimation de la matrice fondamentale à partir des paires de points pa et pb
# en correspondance.
# conversion de kpa et kpb en NumPy array de taille (n,3) contenant les coordonnées (u,v,1) de chaque keypoint
pa = np.array([np.array(kpa[i].pt+(1,)) for i in range(len(kpa))])
pb = np.array([np.array(kpb[i].pt+(1,)) for i in range(len(kpb))])
# On restreint les listes de points pa et pb pour ne garder que les points listés dans "matches"
qa = np.array([pa[m.queryIdx,:] for m in matches])
qb = np.array([pb[m.trainIdx,:] for m in matches])
# estimation de la matrice fondamentale
Fest, _ = cv2.findFundamentalMat(qa,qb,cv2.FM_RANSAC)
print(np.linalg.norm(F-Fest)/np.linalg.norm(F))
%matplotlib inline
plt.figure(figsize=(14,7))
drawEpilines(Fest,qa,Ib)
plt.plot(qb[:,0],qb[:,1],'o')
plt.tight_layout()
plt.axis([0,Ib.shape[1],Ib.shape[0],0])
plt.show()
0.015288171240025746
# Question 8 : On donne ci-dessous les matrices de projection pour Ia et Ib. Calculer et afficher la reconstruction
# 3D des points en correspondance pa et pb.
Ma = 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]])
Mb = 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]])
# Calcul des coordonnées 3D des points
# N.B. la fonction Reconstr3D est donnée au début du ce notebook,
# elle est directement issue du code de la première partie
Xrec = Reconstr3D(qa, qb, Ma, Mb)
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook
fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(111, projection='3d')
ax.plot(Xrec[:,0],Xrec[:,1],Xrec[:,2],'.')
plt.tight_layout()
plt.show()
# Question 4 : Charger à présent les deux images head_0000.ppm et head_0001.ppm,
# ainsi que head_0000_silh.pgm et head_0001_silh.pgm.
# Afficher ces 4 images.
Ia = plt.imread('head_0000.ppm')
Ib = plt.imread('head_0001.ppm')
silha = plt.imread('head_0000_silh.pgm')
silhb = plt.imread('head_0001_silh.pgm')
%matplotlib inline
plt.figure(figsize=(14,14))
plt.subplot(2,2,1)
plt.imshow(Ia)
plt.subplot(2,2,2)
plt.imshow(Ib)
plt.subplot(2,2,3)
plt.imshow(silha, cmap='gray')
plt.subplot(2,2,4)
plt.imshow(silhb, cmap='gray')
plt.tight_layout()
plt.show()
# Question 5 : En combinant les images head_0000 et head_0000_silh, créer une image Ia correspondant à l'image
# de la statue avec un fond noir.
# Créer de même l'image Ib à partir de head_0001 et head_0001_silh
Ia.shape # (768,1024,3)
silha.shape # (768, 1024)
#Ia = Ia * (255-silha[:,:,None])/255
#Ib = Ib * (255-silhb[:,:,None])/255
Ia = Ia * (silha[:,:,None]==0)
Ib = Ib * (silhb[:,:,None]==0)
Ia = np.uint8(Ia)
Ib = np.uint8(Ib)
plt.figure(figsize=(14,7))
plt.subplot(1,2,1)
plt.imshow(Ia)
plt.subplot(1,2,2)
plt.imshow(Ib)
plt.tight_layout()
plt.show()
# Question 6 : Calculer et afficher les points ORB pour Ia et Ib ainsi que les correspondances.
# On notera pa et pb ces deux ensembles de points.
orb = cv2.ORB_create(edgeThreshold = 200)
# Détection des points ORB sur les deux images
kpa = orb.detect(Ia,None)
kpb = orb.detect(Ib,None)
# Affichage des images avec leurs points
resa = cv2.drawKeypoints(Ia,kpa,None,flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
resb = cv2.drawKeypoints(Ib,kpb,None,flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
%matplotlib inline
plt.figure(figsize=(14,7))
plt.subplot(1,2,1)
plt.imshow(resa)
plt.subplot(1,2,2)
plt.imshow(resb)
plt.tight_layout()
plt.show()
# Affichage des images avec leurs points
resa = cv2.drawKeypoints(Ia,kpa,None,flags=0)
resb = cv2.drawKeypoints(Ib,kpb,None,flags=0)
%matplotlib inline
plt.figure(figsize=(14,7))
plt.subplot(1,2,1)
plt.imshow(resa)
plt.subplot(1,2,2)
plt.imshow(resb)
plt.tight_layout()
plt.show()
# Calcul des descripteurs ORB pour les deux images
kpa, dsca = orb.compute(Ia,kpa)
kpb, dscb = orb.compute(Ib,kpb)
# Création d'un objet pour la mise en correspondance des keypoints
#bf = cv2.BFMatcher(cv2.NORM_L2)
bf = cv2.BFMatcher(cv2.NORM_L2)
# Calcul des correspondances
matches = bf.match(dsca, dscb)
# Affichage des correspondances
image_matches = cv2.drawMatches(Ia, kpa, Ib, kpb, matches, None)
plt.figure(figsize=(14,7))
plt.imshow(image_matches)
plt.tight_layout()
plt.show()
# Question 7 : Calculer une estimation de la matrice fondamentale à partir des paires de points pa et pb
# en correspondance.
# conversion de kpa et kpb en NumPy array de taille (n,3) contenant les coordonnées (u,v,1) de chaque keypoint
pa = np.array([np.array(kpa[i].pt+(1,)) for i in range(len(kpa))])
pb = np.array([np.array(kpb[i].pt+(1,)) for i in range(len(kpb))])
# On restreint les listes de points pa et pb pour ne garder que les points listés dans "matches"
qa = np.array([pa[m.queryIdx,:] for m in matches])
qb = np.array([pb[m.trainIdx,:] for m in matches])
# estimation de la matrice fondamentale
F, mask = cv2.findFundamentalMat(qa,qb,cv2.FM_RANSAC)
# Affichage des correspondances
matches = list(np.array(matches)[mask.astype(bool).flatten()])
image_matches = cv2.drawMatches(Ia, kpa, Ib, kpb, matches, None)
%matplotlib inline
plt.figure(figsize=(14,7))
plt.imshow(image_matches)
plt.tight_layout()
plt.show()
plt.figure(figsize=(14,7))
drawEpilines(F,qa,Ib)
plt.plot(qb[:,0],qb[:,1],'o')
plt.axis([0,Ib.shape[1],Ib.shape[0],0])
plt.tight_layout()
plt.show()
# Question 8 : On donne ci-dessous les matrices de projection pour Ia et Ib. Calculer et afficher la reconstruction
# 3D des points en correspondance pa et pb.
Ma = np.array([[-487.257,1099.86,-78.1951,33216.2],
[-179.829,23.2411,-1145.54,22371.4],
[-0.981925,0.0212885,-0.188068,61.5762]])
Mb = np.array([[-732.432,954.848,-70.9875,31937 ],
[-178.143,-9.49754,-1146,22240.6 ],
[-0.958974,-0.208851,-0.191706,61.6066 ]])
# Calcul des coordonnées 3D des points
# N.B. la fonction Reconstr3D est donnée au début du ce notebook,
# elle est directement issue du code de la première partie
Xrec = Reconstr3D(qa, qb, Ma, Mb)
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook
fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(111, projection='3d')
ax.plot(Xrec[:,0],Xrec[:,1],Xrec[:,2],'.')
plt.tight_layout()
plt.show()
len(matches)
72