import numpy as np
import matplotlib.pyplot as plt
# question 1
def plotmesh(X,Y):
# fonction pour afficher une grille de points
# X,Y : matrices de même taille donnant les coordonnées des points de la grille
plt.plot(X,Y,'g');
plt.plot(X.T,Y.T,'g')
# correction
M, N = 50, 25
X, Y = np.meshgrid(np.arange(M),np.arange(N),indexing='ij')
plotmesh(X,Y)
plt.axis('equal')
plt.show()
# question 2
def rotation(X,Y,theta,cx,cy):
# Rotation d'un ensemble de points 2D.
# X,Y : NumPy arrays de même taille donnant les coordonnées des points d'entrée
# theta : angle de la rotation en radians
# cx,cy : coordonnées du centre de la rotation
# Xout,Yout : tableaux de même taille que X et Y donnant les coordonnées des points après rotation.
sz = X.shape
R = np.array([[np.cos(theta),-np.sin(theta)],[np.sin(theta),np.cos(theta)]])
XYout = R @ np.array([X.flatten()-cx,Y.flatten()-cy])
XYout[0,:] += cx
XYout[1,:] += cy
Xout = np.reshape(XYout[0,:],sz)
Yout = np.reshape(XYout[1,:],sz)
return Xout, Yout
# to do : rotation de la grille et affichage
theta = (np.pi/180) * 30
cx, cy = M/2, N/2
RX, RY = rotation(X,Y,theta,cx,cy)
plotmesh(RX,RY)
plt.axis('equal')
plt.show()
# question 3
# to do : projection de la grille tournée sur la grille initiale et affichage
pRX, pRY = np.round(RX), np.round(RY)
pRX, pRY = np.maximum(pRX,0), np.maximum(pRY,0)
pRX, pRY = np.minimum(pRX,M-1), np.minimum(pRY,N-1)
plotmesh(pRX,pRY)
plt.axis('equal')
plt.show()
# question 4
# correction
I1 = plt.imread('Barbara.bmp')
plt.figure(figsize=(10,10))
plt.imshow(I1,cmap='gray')
plt.show()
# question 5
def imrot_v1(I,theta):
M, N = I.shape
J = np.zeros((M,N))
X, Y = np.meshgrid(np.arange(M),np.arange(N),indexing='ij')
cx, cy = M/2, N/2
RX, RY = rotation(X,Y,theta,cx,cy)
pRX, pRY = np.round(RX), np.round(RY)
pRX, pRY = np.maximum(pRX,0), np.maximum(pRY,0)
pRX, pRY = np.minimum(pRX,M-1), np.minimum(pRY,N-1)
pRX, pRY = (pRX).astype(int), (pRY).astype(int)
for i in range(M):
for j in range(N):
J[pRX[i,j],pRY[i,j]] = I[i,j]
return J
theta = np.pi/5
RI1_v1 = imrot_v1(I1,theta)
plt.figure(figsize=(10,10))
plt.imshow(RI1_v1,cmap='gray')
plt.show()
# question 6
def imrot_v2(I,theta):
M, N = I.shape
J = np.zeros((M,N))
X, Y = np.meshgrid(np.arange(M),np.arange(N),indexing='ij')
cx, cy = M/2, N/2
RX, RY = rotation(X,Y,-theta,cx,cy)
pRX, pRY = np.round(RX), np.round(RY)
pRX, pRY = np.maximum(pRX,0), np.maximum(pRY,0)
pRX, pRY = np.minimum(pRX,M-1), np.minimum(pRY,N-1)
pRX, pRY = (pRX).astype(int), (pRY).astype(int)
for i in range(M):
for j in range(N):
J[i,j] = I[pRX[i,j],pRY[i,j]]
return J
RI1_v2 = imrot_v2(I1,theta)
plt.figure(figsize=(10,10))
plt.imshow(RI1_v2,cmap='gray')
plt.show()
# question 7
import scipy.ndimage
def imrot(I, theta, order=1):
M, N = I.shape
X, Y = np.meshgrid(np.arange(M),np.arange(N),indexing='ij')
cx, cy = M/2, N/2
RX, RY = rotation(X,Y,-theta,cx,cy)
coordinates = [RX.flatten(), RY.flatten()]
J = scipy.ndimage.map_coordinates(I, coordinates, order=order, mode='constant')
J = J.reshape((M,N))
return J
RI1 = imrot(I1,theta)
plt.figure(figsize=(10,10))
plt.imshow(RI1,cmap='gray')
plt.show()
# version avec interpolation bicubique au lieu de bilinéaire.
# On ne remarque pas de différences sur cet exemple.
RI1 = imrot(I1,theta, order=3)
plt.figure(figsize=(10,10))
plt.imshow(RI1,cmap='gray')
plt.show()
# question 8
theta = np.pi/6
RI1 = I1.copy()
for i in range(12):
RI1 = imrot(RI1,theta)
plt.figure(figsize=(10,10))
plt.imshow(RI1,cmap='gray')
plt.show()
theta = np.pi/6
RI1 = I1.copy()
for i in range(12):
RI1 = imrot(RI1,theta, order=3)
plt.figure(figsize=(10,10))
plt.imshow(RI1,cmap='gray')
plt.show()
# question 9
def homothetie(X,Y,rho,cx,cy):
# Homothétie d'un ensemble de points 2D.
# X,Y : NumPy arrays de même taille donnant les coordonnées des points d'entrée
# rho : facteur de l'homothétie
# cx,cy : coordonnées du centre de la rotation
# Xout,Yout : tableaux de même taille que X et Y donnant les coordonnées des points après rotation.
Xout = rho*(X-cx) + cx
Yout = rho*(Y-cy) + cy
return Xout, Yout
def imzoom(I, rho, order=1):
M, N = I.shape
X, Y = np.meshgrid(np.arange(M),np.arange(N),indexing='ij')
cx, cy = M/2, N/2
RX, RY = homothetie(X,Y,1/rho,cx,cy)
coordinates = [RX.flatten(), RY.flatten()]
J = scipy.ndimage.map_coordinates(I, coordinates, order=order, mode='constant')
J = J.reshape((M,N))
return J
rho = 1.5
RI1 = imzoom(I1,rho)
plt.figure(figsize=(10,10))
plt.imshow(RI1,cmap='gray')
plt.show()
rho = 0.5
RI1 = imzoom(I1,rho)
plt.figure(figsize=(10,10))
plt.imshow(RI1,cmap='gray')
plt.show()
def transfo_affine(X,Y, M, tau):
# transofrmation affine d'un ensemble de points 2D.
# X,Y : NumPy arrays de même taille donnant les coordonnées des points d'entrée
# M, tau : matrie et vecteur de translation
# Xout,Yout : tableaux de même taille que X et Y donnant les coordonnées des points après rotation.
sz = X.shape
XYout = M @ np.array([X.flatten(),Y.flatten()])
XYout[0,:] += tau[0]
XYout[1,:] += tau[1]
Xout = np.reshape(XYout[0,:],sz)
Yout = np.reshape(XYout[1,:],sz)
return Xout, Yout
def im_transfo_affine(I, A, tau, order=1, inc_out = 1):
M, N = I.shape
Mout, Nout = M*inc_out, N*inc_out
X, Y = np.meshgrid(np.arange(Mout),np.arange(Nout),indexing='ij')
Ainv = np.linalg.inv(A)
RX, RY = transfo_affine(X,Y, Ainv,-Ainv@tau)
coordinates = [RX.flatten(), RY.flatten()]
J = scipy.ndimage.map_coordinates(I, coordinates, order=order, mode='constant')
J = J.reshape((Mout,Nout))
return J
A = np.array([[0.4,0.2],[-0.2,0.9]])
tau = np.array([[50],[-10]])
RI1 = im_transfo_affine(I1, A, tau)
plt.figure(figsize=(10,10))
plt.imshow(RI1,cmap='gray')
plt.show()
import scipy.io
tmp = scipy.io.loadmat('field.mat')
tmp['Uinv'].shape
(512, 512)
# question 12
import scipy.io
def imdef_field(I,Uinv,Vinv, order=1):
M, N = I.shape
X, Y = np.meshgrid(np.arange(M),np.arange(N),indexing='ij')
Xinv = X + Uinv
Yinv = Y + Vinv
coordinates = [Xinv.flatten(), Yinv.flatten()]
J = scipy.ndimage.map_coordinates(I, coordinates, order=order, mode='constant')
J = J.reshape((M,N))
return J
tmp = scipy.io.loadmat('field.mat')
Uinv, Vinv = tmp['Uinv'], tmp['Vinv']
I1d = imdef_field(I1,Uinv,Vinv, order=3)
plt.figure(figsize=(10,10))
plt.imshow(I1d,cmap='gray')
plt.show()
# question 13
def invchamp(Uinv,Vinv, order=1):
M, N = Uinv.shape
X, Y = np.meshgrid(np.arange(M),np.arange(N),indexing='ij') # correspond à y dans la formule d'itérations
niter = 10
Xi, Yi = X.copy(), Y.copy() # correspond à x^n dans la formule
for i in range(niter):
# on doit calculer le champ Uinv, Vinv interpolé sur la grille déformée Xi, Yi
coordinates = [Xi.flatten(), Yi.flatten()]
W1 = scipy.ndimage.map_coordinates(Uinv, coordinates, order=order, mode='constant')
W1 = W1.reshape((M,N))
Xi = X - W1
W2 = scipy.ndimage.map_coordinates(Vinv, coordinates, order=order, mode='constant')
W2 = W2.reshape((M,N))
Yi = Y - W2
# à ce stade Xi, Yi correpondent au point x^* limite.
U = X-Xi
V = Y-Yi
return U, V
U, V = invchamp(Uinv,Vinv)
I1_ = imdef_field(I1d,U,V)
plt.figure(figsize=(10,10))
plt.imshow(I1_,cmap='gray')
plt.show()