# script pour le recalage d'images basé sur des points de référence
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [10, 10]
import scipy.ndimage
def PlotPoints(x,clr='b',num=None):
# affichage d'un liste de points numerotes
# x : matrice 2*n des coordonnees des points
# clr : optionnel, caractere donnant la couleur voulue : 'r', 'b', etc.
# par defaut clr = 'b'
# num : optionnel, vecteur d'entiers contenant les numeros
# par defaut num = 0:n-1
n = x.shape[1]
if num==None:
num = range(n)
for k in range(n):
# affichage du point
plt.plot(x[0,k],x[1,k],'x'+clr)
# affichage du numero du point
plt.text(x[0,k]+3,x[1,k]-10,num[k])
def affine_trans(X,Y,A,b):
# transformation affine d'un ensemble de points 2D
# entrées:
# X,Y : tableaux de meme taille donnant les coordonnees des points d'entree
# A : matrice 2*2 de la partie lineaire de la transformation
# b : vecteur 2*1 de translation
# sorties:
# Xout,Yout : tableaux de meme taille que X et Y donnant les coordonnees
# des points apres transformation.
sz = X.shape
XYout = A @ np.array([X.flatten(),Y.flatten()])
XYout[0,:] += b[0]
XYout[1,:] += b[1]
Xout = np.reshape(XYout[0,:],sz)
Yout = np.reshape(XYout[1,:],sz)
return Xout, Yout
def imdef_affine(I,A,b):
# 0. définition de la grille de pixels
m = I.shape[1]
n = I.shape[0]
X, Y = np.meshgrid(range(m),range(n))
# 1. calcul des images des points de la grille par la transformation affine inverse
Ainv = np.linalg.inv(A)
Xnew, Ynew = affine_trans(X,Y,Ainv,-Ainv@b)
# 2. conversion de Xnew et Ynew en tableau P à 2 lignes
# P doit être de taille 2*N, et donner les coordonnées des pixels
# de la grille ayant subi la rotation
P = np.array([Ynew.flatten(),Xnew.flatten()])
# 3. interpolation
Inew = scipy.ndimage.map_coordinates(I,P,order=1)
# ici Inew est de taille N ; il faut donc changer sa taille à nouveau
Inew = Inew.reshape(I.shape)
return Inew
# chargement des données
import scipy.io as sio
fish = sio.loadmat('fish.mat') # contient les images I et J et leurs points x et y
I = fish['I']
J = fish['J']
x = fish['x']
y = fish['y']
# On affiche les deux images et leurs points
plt.rcParams['figure.figsize'] = [16, 8]
plt.subplot(1,2,1)
plt.imshow(I, cmap=plt.cm.gray)
PlotPoints(x)
plt.subplot(1,2,2)
plt.imshow(J, cmap=plt.cm.gray)
PlotPoints(y)
plt.show()
fish
{'__header__': b'MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Mon Nov 13 00:53:05 2017', '__version__': '1.0', '__globals__': [], 'I': array([[255., 255., 255., ..., 255., 255., 255.], [255., 255., 255., ..., 255., 255., 255.], [255., 255., 255., ..., 255., 255., 255.], ..., [255., 255., 255., ..., 255., 255., 255.], [255., 255., 255., ..., 255., 255., 255.], [255., 255., 255., ..., 255., 255., 255.]]), 'J': array([[255., 255., 255., ..., 255., 255., 255.], [255., 255., 255., ..., 255., 255., 255.], [255., 255., 255., ..., 255., 255., 255.], ..., [255., 255., 255., ..., 255., 255., 255.], [255., 255., 255., ..., 255., 255., 255.], [255., 255., 255., ..., 255., 255., 255.]]), 'x': array([[ 46.35267233, 132.23946616, 201.36362338, 254.9330723 , 317.89884668, 388.74416147, 352.39649552, 348.32147356, 298.47725953, 234.13523959, 164.14806605, 93.17048928, 34.65996781, 76.09936409, 71.81242843, 120.35213926, 126.1693807 , 87.50804437, 66.01661706, 114.65363909, 156.51309217, 124.12060245, 110.40202238, 197.97539115], [ 31.94817912, 49.71683054, 83.32309968, 118.32892158, 168.59749813, 183.96423808, 207.59226335, 250.2552956 , 195.64374559, 176.03540168, 169.28790093, 117.9613759 , 43.17904477, 40.24107394, 48.53266085, 58.39580146, 80.7005735 , 99.00266847, 80.4799588 , 91.61904678, 106.03618373, 118.33894986, 105.63920791, 41.52922893]]), 'y': array([[210.76136417, 203.17941908, 187.9744165 , 246.86535298, 317.18538616, 337.38287723, 358.98009875, 371.39124819, 340.26876878, 362.17784666, 365.81919978, 307.69610588, 218.53275738, 221.0654437 , 238.64101392, 238.68352826, 258.08724925, 267.20884885, 251.94498156, 262.13226065, 268.43244118, 282.20757997, 275.99242812, 156.66085055], [405.99927045, 350.28198376, 280.61392863, 265.95864337, 274.2637878 , 231.65095356, 240.75901939, 257.9698785 , 289.79817027, 329.26293805, 399.30091004, 434.52454757, 418.04344035, 381.56271108, 395.23288748, 366.21252025, 378.96775104, 408.02243309, 426.4964574 , 379.01468339, 312.76467341, 368.40979268, 376.6190682 , 237.33142389]])}
# calcul des transformations optimales
# methode 0 : recalage par translation uniquement
def FindOptimTrans(x,y):
# calcul de la translation optimale
# entre deux ensembles de points xi et yi dans R^d
# entrées:
# x, y : tailles d*n, coordonnees des points
# sorties :
# A : matrice identite d*d, matrice de la partie lineaire
# b : vecteur d*1, vecteur translation
d = x.shape[0]
A = np.identity(d)
b = np.mean(y,1) - np.mean(x,1)
return A, b
A, b = FindOptimTrans(x,y)
# calcul de l'image transformée
Inew = imdef_affine(I,A,b)
# affichage de l'image transformée
plt.imshow(Inew, cmap=plt.cm.gray)
PlotPoints(A@x+b.reshape((2,1)))
plt.show()
# On affiche une image en couleurs contenant l'image
# et les points transformes en bleu, et l'image et les points cible en rouge
K = np.stack((Inew/256,J/256,np.ones(I.shape)),axis=2)
plt.imshow(K)
PlotPoints(y,'r')
z1, z2 = affine_trans(x[0,:],x[1,:],A,b)
z = np.vstack((z1,z2))
PlotPoints(z,'b')
plt.plot(np.vstack((y[0,:],z[0,:])),np.vstack((y[1,:],z[1,:])),'b')
plt.show()
np.hstack((np.random.rand(2,2),np.random.rand(2,2)))
array([[0.96577089, 0.08152954, 0.61176285, 0.26916205], [0.06559388, 0.66652328, 0.69771119, 0.52421374]])
# methode 1 : recalage par transformations affines
def FindOptimAffine(x,y):
# calcul de la transfo affine optimale
# entre deux ensembles de points xi et yi dans R^d
# entrées:
# x, y : tailles d*n, coordonnees des points
# sorties :
# A : matrice d*d, matrice de la partie lineaire
# b : vecteur d*1, vecteur translation
d, n = x.shape
xtilde = np.vstack((x,np.ones((1,n))))
Mt = np.linalg.solve(xtilde@xtilde.T, xtilde@y.T) # ou bien : np.linalg.inv(xtilde@xtilde.T) @ xtilde@y.T
M = Mt.T
A = M[:,:d]
b = M[:,d]
return A, b
A, b = FindOptimAffine(x,y)
# calcul de l'image transformée
Inew = imdef_affine(I,A,b)
# affichage de l'image transformée
plt.imshow(Inew, cmap=plt.cm.gray)
plt.show()
# On affiche une image en couleurs contenant l'image
# et les points transformes en bleu, et l'image et les points cible en rouge
K = np.stack((Inew/256,J/256,np.ones(I.shape)),axis=2)
plt.imshow(K)
PlotPoints(y,'r')
z1, z2 = affine_trans(x[0,:],x[1,:],A,b)
z = np.vstack((z1,z2))
PlotPoints(z,'b')
plt.plot(np.vstack((y[0,:],z[0,:])),np.vstack((y[1,:],z[1,:])),'b')
plt.show()
# methode 2 : recalage par déplacements
def FindOptimDeplace(x,y):
# calcul du déplacement optimal
# entre deux ensembles de points xi et yi dans R^2
# entrées:
# x, y : tailles 2*n, coordonnees des points
# sorties :
# A : matrice de la partie lineaire
# b : vecteur 2*1, vecteur translation
# ... to do ...
return A, b
A, b = FindOptimDeplace(x,y)
# calcul de l'image transformée
Inew = imdef_affine(I,A,b)
# affichage de l'image transformée
plt.imshow(Inew, cmap=plt.cm.gray)
plt.show()
# On affiche une image en couleurs contenant l'image
# et les points transformes en bleu, et l'image et les points cible en rouge
K = np.stack((Inew/256,J/256,np.ones(I.shape)),axis=2)
plt.imshow(K)
PlotPoints(y,'r')
z1, z2 = affine_trans(x[0,:],x[1,:],A,b)
z = np.vstack((z1,z2))
PlotPoints(z,'b')
plt.plot(np.vstack((y[0,:],z[0,:])),np.vstack((y[1,:],z[1,:])),'b')
plt.show()
# methode 3 : homographies
# méthode 4 : Thin Plate Splines
def FindOptimTPS(x,y):
# calcul de la transformation TPS optimale
# entre deux ensembles de points xi et yi dans R^2
# entrées:
# x, y : tailles 2*n, coordonnées des points
# sorties :
# alpha : taille 2*x, coordonnées des vecteurs alpha_i de la transformation TPS
# A : matrice de la partie linéaire
# b : vecteur 2*1, vecteur translation
d, n = x.shape # d=2 normalement$
ones = np.ones((1,n))
xtilde = np.vstack((x,ones))
# calcul de la matrice h(x) contenant les h(||x_i-x_j||)
def h(r):
return r**2 * np.log(r)
# calcul des x_i - x_j : x (2,n)->(2,n,1) et x (2,n)->(2,1,n)
xmx = x[:,:,None] - x[:,None,:]
# calcul des ||x_i-x_j||
norm_xmx = np.sqrt(np.sum(xmx**2,axis=0))
# calcul des h(||x_i-x_j||)
hx = h(norm_xmx)
return alpha, A, b