import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [10, 10]
import scipy.ndimage
/Users/glaunes/opt/miniconda3/lib/python3.9/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.1 warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
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
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()
A = np.array([[1.,-.1],[.2,1.]])
b = np.array([[20.],[20.]])
# 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()
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 :
# tau : vecteur d*1, vecteur translation
tau = np.mean(y,axis=1) - np.mean(x,axis=1)
return tau
tau = FindOptimTrans(x,y)
A, b = np.identity(2), tau
# 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()
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
# ... to do ...
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()
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()
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
# ... to do ...
return alpha, A, b