# 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()]) + b
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
n, m = I.shape
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()
# 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))[:,None]
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()
# 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 identite d*d, matrice de la partie lineaire
# b : vecteur d*1, vecteur translation
d, n = x.shape
ones = np.ones((1,n))
xtilde = np.vstack((x,ones)) # ou np.concatenate((x,ones),axis=0)
ytilde = np.vstack((y,ones)) # ou np.concatenate((y,ones),axis=0)
mat_x = xtilde @ xtilde.T
mat_xy = y @ xtilde.T
# 1e solution : utilisation de np.linalg.solve pour résoudre (xtilde@xtilde.T).T @ M.T = (y@xtilde.T).T
# M = np.linalg.solve(mat_x.T,mat_xy.T).T
# 2e solution : M @ xtilde@xtilde.T = y@xtilde.T => M = y@xtilde.T @ (xtilde@xtilde.T)^-1
M = mat_xy @ np.linalg.inv(mat_x)
A = M[:,:d]
b = M[:,d][:,None]
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)
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()
# 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
n = x.shape[1]
x = x[0,:] + 1j * x[1,:]
y = y[0,:] + 1j * y[1,:]
sx = np.sum(x)
sy = np.sum(y)
syx = np.sum(y*np.conjugate(x))
alpha = sy*np.conjugate(sx)/n - syx
r = np.abs(alpha)
nu = np.angle(alpha)
if r==0:
theta = 0
b = (sy - sx)/n
else:
theta1 = nu
b1 = (sy - sx * np.exp(1j*theta1))/n
J1 = np.sum(np.abs(np.exp(1j*theta1)*x+b1-y)**2)
theta2 = nu + np.pi
b2 = (sy - sx * np.exp(1j*theta2))/n
J2 = np.sum(np.abs(np.exp(1j*theta2)*x+b2-y)**2)
if J1 < J2:
theta = theta1
b = b1
else:
theta = theta2
b = b2
A = np.array([[np.cos(theta),-np.sin(theta)],
[np.sin(theta),np.cos(theta)]])
b = np.array([[np.real(b)],[np.imag(b)]])
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)
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()
# methode 3 : homographies
# méthode 4 : Thin Plate Splines
def x2logx(r):
# Fonction h(r) = r^2 log(r) (avec h(0)=0) vectorisée.
# r est un numpy.array de taille arbitraire
# N.B. On évite de calculer des log(0) à l'aide d'un masque binaire sélectionnant uniquement
# les valeurs non nulles dans r.
res = r**2
mask = r>0
res[mask] *= np.log(r[mask])
return res
def KernelMatrix(h, x, y):
# Calcul d'une matrice noyau.
# entrées : h fonction scalaire vectorisée
# x tableau NumPy de taille (d,m), coordonnées de m points dans R^d
# y tableau NumPy de taille (d,n), coordonnées de n points dans R^d
# sortie : H tableau NumPy de taille (m,n) tel que H_ij = h(|x_i-y_j|) où |.| est le norme euclidienne
#
# calcul des x_i - y_j :
xmy = x[:,:,None] - y[:,None,:] # N.B. xmy est de taille (d,m,n)
# calcul des ||x_i-y_j||
norm_xmy = np.sqrt(np.sum(xmy**2,axis=0))
# calcul des h(||x_i-y_j||)
H = h(norm_xmy)
return H
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,n), coordonnées des vecteurs alpha_i de la transformation TPS
# A : taille (2,2), matrice de la partie linéaire
# b : taille (2,1), vecteur translation
d, n = x.shape # d=2 normalement
ones = np.ones((1,n))
xtilde = np.vstack((x,ones))
# On calcule de la matrice H de taille (n,n) contenant les h(||x_i-x_j||)
H = KernelMatrix(x2logx,x,x)
# On forme la matrice de taille (n+3,n+3) du système linéaire à résoudre (cf document pdf)
K = np.block([[H,xtilde.T],[xtilde,np.zeros((3,3))]])
# On forme le tableau de taille (n+3,2), second membre du système linéaire
Y = np.vstack((y.T,np.zeros((3,2))))
# On résoud le système linéaire
res = np.linalg.solve(K,Y)
# On extrait alpha, A et b du résultat
alpha = res[:n,:].T
A = res[n:n+2,:].T
b = res[n+2,:][:,None] # N.B : res[n+2,:] est de taille (2), res[n+2,:][:,None] est de taille (2,1)
return alpha, A, b
def TPS_trans(X,Y,x,alpha,A,b):
# transformation Thin Plate Splines (TPS) d'un ensemble de points 2D
# entrées:
# X,Y : tableaux de meme taille donnant les coordonnees des points d'entrée
# x : tableau de taille (2,n), points de référence de la transformation TPS,
# alpha : tableau de taille (2,n), coefficients de la transformation TPS,
# A : matrice 2*2 de la partie lineaire de la transformation TPS,
# b : vecteur 2*1 de translation de la transformation TPS,
# sorties:
# Xout,Yout : tableaux de meme taille que X et Y donnant les coordonnees
# des points apres transformation.
sz = X.shape
XY = np.array([X.flatten(),Y.flatten()])
H = KernelMatrix(x2logx,XY,x)
XYout = (H @ alpha.T).T + A @ XY + b
Xout = np.reshape(XYout[0,:],sz)
Yout = np.reshape(XYout[1,:],sz)
return Xout, Yout
def imdef_TPS(I,x,alpha,A,b):
# 0. définition de la grille de pixels
n, m = I.shape
X, Y = np.meshgrid(range(m),range(n))
# 1. calcul des images des points de la grille par la transformation TPS
Xnew, Ynew = TPS_trans(X,Y,x,alpha, A, 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 transformation
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
alpha, A, b = FindOptimTPS(y,x)
# calcul de l'image transformée
Inew = imdef_TPS(I,y,alpha,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
# transformée en rouge, et l'image et les points cible en bleu
K = np.stack((Inew/256,J/256,np.ones(I.shape)),axis=2)
plt.imshow(K)
PlotPoints(y,'r')
plt.show()