import numpy as np
import matplotlib.pyplot as plt
# Question 2
h1 = np.zeros((512,512))
h1[246:346,150:350] = 1
plt.imshow(h1,cmap='gray',vmin=0,vmax=1)
plt.show()
# question 3
# version vue pendant le cours
def ImTranslate(h,tau):
M, N = h.shape
res = np.zeros((M,N))
res[tau[0]:M,tau[1]:N] = h[0:(M-tau[0]),0:(N-tau[1])]
res[0:tau[0],0:tau[1]] = h[(M-tau[0]):M,(N-tau[1]):N]
res[0:tau[0],tau[1]:N] = h[(M-tau[0]):M,0:(N-tau[1])]
res[tau[0]:M,0:tau[1]] = h[0:(M-tau[0]),(N-tau[1]):N]
return res
# version alternative (avec l'utilisation de l'opération modulo)
def ImTranslate(h,tau):
M, N = h.shape
res = np.zeros((M,N))
indrow = ((tau[0]+np.arange(M))%M).reshape(M,1)
indcol = ((tau[1]+np.arange(N))%N).reshape(1,N)
res[indrow,indcol] = h
return res
h1_tr = ImTranslate(h1,[250,250])
plt.imshow(h1_tr,cmap='gray',vmin=0,vmax=1)
plt.show()
# question 4
def DSN(h,n):
M, N = h.shape
res = 0
for i in range(n):
tau = np.random.randint(M), np.random.randint(N)
res += ImTranslate(h,tau)
return res
# on fait une fonction pour afficher l'image et ses DSN avec des valeurs croissantes de n:
def ShowDSN(h):
ns = [10, 100, 1000]
plt.figure(figsize=(16,8))
plt.subplot(1,len(ns)+1,1)
plt.imshow(h,cmap='gray')
plt.title('original image')
for i, n in enumerate(ns):
h_DSN = DSN(h,n)
plt.subplot(1,len(ns)+1,i+2)
plt.imshow(h_DSN,cmap='gray')
plt.title('DSN, n='+str(n))
plt.show()
ShowDSN(h1)
h2 = plt.imread('images/coeur.bmp')
h2 = np.mean(h2,2)
ShowDSN(h2)
h3 = plt.imread('images/briques.bmp')
h3 = np.mean(h3,2)
ShowDSN(h3)
h4 = plt.imread('images/mer.bmp')
h4 = np.mean(h4,2)
ShowDSN(h4)
# question 5
def GaussImage(N,b=20):
if isinstance(b,int) or isinstance(b,float):
b *= np.array([-1,1,-1,1])
h = np.exp(-(np.ones((N,1))*np.linspace(b[0],b[1],N)[None,:])**2
-(np.ones((N,1))*np.linspace(b[2],b[3],N)[None,:]).T**2)
return h
def ShowDFT(h):
# DFT de l'image h
fh = np.fft.fft2(h)
# module de la DFT
afh = np.abs(fh)
# phase de la DFT
pfh = np.angle(fh)
# affichage
plt.figure(figsize=(16,8))
plt.subplot(1,3,1)
plt.imshow(h,cmap='gray')
plt.title('image')
plt.subplot(1,3,2)
plt.imshow(np.log(0.1+np.fft.fftshift(afh)),cmap="gray")
plt.title('amplitude de la DFT (centrée et échelle log)')
plt.subplot(1,3,3)
plt.imshow(np.fft.fftshift(pfh),cmap="gray")
plt.title('phase de la DFT (centrée)')
plt.show()
# l'image du rectangle
ShowDFT(h1)
# quelques images du dossier
ShowDFT(h2)
ShowDFT(h3)
ShowDFT(h4)
# image d'une gaussienne
h5 = GaussImage(512,40)
ShowDFT(h5)
# image de bruit blanc
h6 = np.random.randn(512,512)
ShowDFT(h6)
# question 6
def ADSN(h):
M, N = h.shape
h_tilde = (h-np.mean(h)) / np.sqrt(M*N)
X = np.random.randn(M,N)
Y = np.fft.ifft2( np.fft.fft2(h_tilde) * np.fft.fft2(X) )
Y = np.real(Y)
return Y
# on fait une fonction pour afficher l'image, son DSN avec n grand, et son ADSN
def ShowADSN(h):
ns = [10, 100, 1000]
plt.figure(figsize=(16,8))
plt.subplot(1,3,1)
plt.imshow(h,cmap='gray')
plt.title('original image')
plt.subplot(1,3,2)
plt.imshow(DSN(h,1000),cmap='gray')
plt.title('DSN, n=1000')
plt.subplot(1,3,3)
plt.imshow(ADSN(h),cmap='gray')
plt.title('ADSN')
plt.show()
ShowADSN(h1)
ShowADSN(h2)
ShowADSN(h3)
ShowADSN(h4)
# question 7
Y5 = ADSN(h5)
ShowDFT(h5)
ShowDFT(Y5)
# question 8
def RPN(h):
M, N = h.shape
h_tilde = (h-np.mean(h)) / np.sqrt(M*N)
X = np.random.randn(M,N)
theta_fX = np.angle(np.fft.fft2(X))
Z = np.fft.ifft2( np.fft.fft2(h_tilde) * np.exp(1j*theta_fX) )
Z = np.real(Z)
return Z
# on fait une fonction pour afficher l'image, son ADSN et son RPN
def ShowRPN(h):
ns = [10, 100, 1000]
plt.figure(figsize=(16,8))
plt.subplot(1,3,1)
plt.imshow(h,cmap='gray')
plt.title('original image')
plt.subplot(1,3,2)
plt.imshow(ADSN(h),cmap='gray')
plt.title('ADSN')
plt.subplot(1,3,3)
plt.imshow(RPN(h),cmap='gray')
plt.title('RPN')
plt.show()
ShowRPN(h1)
ShowRPN(h2)
ShowRPN(h3)
ShowRPN(h4)
# question 9
def synthese_texture_etendue(h, M_out, N_out, methode_synthese=RPN):
M, N = h.shape
m = np.mean(h)
h_et = np.zeros((M_out,N_out)) + m
h_et[:M,:N] = h
h_et = np.fft.fftshift(h_et)
return methode_synthese(h_et)
for h in [h1,h2,h3,h4]:
plt.figure(figsize=(16,8))
plt.imshow(synthese_texture_etendue(h,1024,1024), cmap='gray')
plt.show()
h7 = plt.imread('images/denim.bmp')
h7 = np.mean(h7,2)
ShowRPN(h7)
plt.figure(figsize=(16,8))
plt.imshow(synthese_texture_etendue(h7,1024,1024), cmap='gray')
plt.show()
h8 = plt.imread('images/cheveux.bmp')
h8 = np.mean(h8,2)
ShowRPN(h8)
plt.figure(figsize=(16,8))
plt.imshow(synthese_texture_etendue(h8,1024,1024), cmap='gray')
plt.show()
# question 10
def ADSN_couleur(h):
M, N, _ = h.shape
fX = np.fft.fft2(np.random.randn(M,N))
Y = np.zeros((M,N,3))
for k in range(3):
hk = h[:,:,k]
hk_tilde = (hk-np.mean(hk)) / np.sqrt(M*N)
Yk = np.fft.ifft2( np.fft.fft2(hk_tilde) * fX )
Y[:,:,k] = np.real(Yk)
Y[:,:,k] += np.mean(h[:,:,k])
Y = np.minimum(1,np.maximum(0,Y))
return Y
def ShowADSN_couleur(h):
ns = [10, 100, 1000]
plt.figure(figsize=(16,8))
plt.subplot(1,2,1)
plt.imshow(h)
plt.title('original image')
plt.subplot(1,2,2)
sh = ADSN_couleur(h)
plt.imshow(sh)
plt.title('ADSN')
plt.show()
u1 = plt.imread('images/cheveux.bmp') / 255
u2 = plt.imread('images/briques.bmp') / 255
u3 = plt.imread('images/mer.bmp') / 255
u4 = plt.imread('images/herbe.bmp') / 255
for u in [u1, u2, u3, u4]:
ShowADSN_couleur(u)
# question 10 avec RPN
def RPN_couleur(h):
M, N, _ = h.shape
theta_fX = np.angle(np.fft.fft2(np.random.randn(M,N)))
Z = np.zeros((M,N,3))
for k in range(3):
hk = h[:,:,k]
hk_tilde = (hk-np.mean(hk)) / np.sqrt(M*N)
Zk = np.fft.ifft2( np.fft.fft2(hk_tilde) * np.exp(1j*theta_fX) )
Z[:,:,k] = np.real(Zk)
Z[:,:,k] *= np.sqrt(M*N)
Z[:,:,k] += np.mean(h[:,:,k])
Z = np.minimum(1,np.maximum(0,Z))
return Z
def ShowRPN_couleur(h):
ns = [10, 100, 1000]
plt.figure(figsize=(16,8))
plt.subplot(1,2,1)
plt.imshow(h)
plt.title('original image')
plt.subplot(1,2,2)
sh = RPN_couleur(h)
plt.imshow(sh)
plt.title('RPN')
plt.show()
for u in [u1, u2, u3, u4]:
ShowRPN_couleur(u)
# Décomposition d'une image en composantes périodique et régulière
def perdecomp(u):
M, N = u.shape
v1 = np.zeros((M,N))
diff1 = u[-1,:] - u[0,:]
v1[0,:] = diff1
v1[-1,:] = - diff1
v2 = np.zeros((M,N))
diff2 = u[:,-1] - u[:,0]
v2[:,0] = diff2
v2[:,-1] = - diff2
v = v1 + v2
fv = np.fft.fft2(v)
X = np.arange(M)
Y = np.arange(N)
fx = np.cos(2*np.pi*X[:,None]/M)
fy = np.cos(2*np.pi*Y[None,:]/N)
div = 2*fx+2*fy-4
div[0,0] = 1
fu = np.fft.fft2(u)
fp = fu - fv / div
fp[0,0] = fu[0,0]
p = np.real(np.fft.ifft2(fp))
s = u - p
return p, s
h10 = plt.imread('images/Barbara.bmp')
h10 = h10 / 255.
h8_01 = h8 / 255
for h in [h10, h8_01]:
p, s = perdecomp(h)
plt.figure(figsize=(16,8))
plt.subplot(1,3,1)
plt.imshow(h,cmap='gray', vmin=0, vmax=1)
plt.subplot(1,3,2)
plt.imshow(p,cmap='gray', vmin=0, vmax=1)
plt.subplot(1,3,3)
plt.imshow(s,cmap='gray')
plt.figure(figsize=(16,8))
plt.subplot(1,2,1)
plt.imshow(np.tile(h,(2,2)),cmap='gray', vmin=0, vmax=1)
plt.subplot(1,2,2)
plt.imshow(np.tile(p,(2,2)),cmap='gray', vmin=0, vmax=1)
plt.show()
# Application à la stnthèse de textures
def perdecomp_couleur(u):
p = np.zeros(u.shape)
s = np.zeros(u.shape)
for k in range(u.shape[2]):
p[:,:,k], s[:,:,k] = perdecomp(u[:,:,k])
return p, s
for u in [u1, u2, u3, u4]:
p, s = perdecomp_couleur(u)
p = np.minimum(1,np.maximum(0,p))
ShowRPN_couleur(p)
u5 = plt.imread('images/blue.bmp') / 255
u6 = plt.imread('images/bois.bmp') / 255
u7 = plt.imread('images/cubes.bmp') / 255
u8 = plt.imread('images/couleurs.bmp') / 255
u9 = plt.imread('images/zebra.bmp') / 255
for u in [u5, u6, u7, u8, u9]:
p, s = perdecomp_couleur(u)
p = np.minimum(1,np.maximum(0,p))
ShowRPN_couleur(p)