Processing math: 100%

TP4 - Suite

On récupère les données de la première partie.

import matplotlib.pyplot as plt
import numpy as np

import matplotlib.image as mpimg

# image
name = 'hibiscus256.bmp'
img0 = mpimg.imread(name)
n,N,a = img0.shape # n : nb de lignes et N : nb de colonnes

# signal
Sr = img0[50,:,0]/255
Sv = img0[50,:,1]/255
Sb = img0[50,:,2]/255
S = (Sr + Sv + Sb)/3

# variations 1D
def delta(V):
    return V[1:] - V[:-1]

Exercice 4.- Descente de gradient pour le problème de débruitage 1D

On a à présent tous les outils permettant de résoudre numériquement le problème $(P)$ dans le cas $q =2$. $$\tag{$P$} \min_{V \in \mathbb{R}^N} E_q(V) \quad \text{avec} \quad E_q(V) = \frac{1}{2} \| V - S_{noise} \|_2^2 + \frac{\lambda}{2} J_q (V) \quad \text{pour} \quad \lambda > 0 \text{ fixé.}$$ On reprend nos signaux $S$ et $S_{noise}$ initiaux et on rappelle qu'on prendra $S_0 = S_{noise}$ dans le terme d'attache aux données de $(P)$.

noise = np.random.normal(0, 0.1, N)
Snoise = S + noise

1. Définir une fonction gradE(V, lam) qui prend en argument un np.array V et le paramètre lam ($\lambda$ dans le modèle) et renvoie $\nabla E(V)$.

En reprenant le résultat obtenu à l'exercice $2$, on remarque que $$\nabla E(V) = V - S_0 + \frac{\lambda}{2}\nabla J(V) \quad \text{et} \quad \frac{1}{2}\nabla J(V) = \left[ V_0 - V_1 \right. ,\underbrace{ - \delta(\delta V) }_{\in \mathbb{R}^{N-2}}, \left. V_{N-1} - V_{N-2} \right] \: ,$$ où $\delta V$ a été défini à l'exercice $1$.

# gradient de E avec parametre lam
def gradE(V,lam):
    N = V.size
    W = np.zeros(N)
    W[1:-1] = - delta(delta(V))
    W[0] = V[0] - V[1]
    W[N-1] = V[N-1] - V[N-2]
    return V - Snoise + lam*W

1. Modifier la fonction gradDescent de l'exercice $3$ en une fonction gradDescentParam(gradfun, lam, x0, nbIter, tau) afin de prendre en compte le paramètre lam apparaissant dans $E$ et dans son gradient gradE(V, lam).

# descente de gradient avec parametre lam
def gradDescentParam(gradfun, lam, x0, nbIter, tau):
    xk = x0
    for k in range(nbIter):
        xk = xk - tau*gradfun(xk, lam)
    return xk

1. On a établi dans l'exercice $2$ que la Hessienne de $E$ était $\displaystyle \nabla^2 E(V) = I_N + \lambda A$. On admettra que $A$ est une matrice symétrique positive dont les valeurs propres sont toutes dans $[0, 4]$. Pourquoi peut-on choisir un pas de temps $\displaystyle 0 < \tau < \frac{2}{1 + 4 \lambda}$ ?

Les valeurs propres de la Hessienne $\displaystyle \nabla^2 E(V) = I_N + \lambda A$ sont exactement $$\{ 1 + \lambda \mu \: : \: \mu \in {\rm Sp}(A) \} \: ,$$ On est dans les hypothèses du théorème de convergence de la descente de gradient à pas fixe des notes de cours avec $l = 1$ et $L = 1 + 4\lambda$, à condition de choisir un pas $\displaystyle 0 < \tau < \frac{2}{L} = \frac{2}{1 + 4 \lambda}$.

On prendra par la suite un pas de temps $\displaystyle \tau = \frac{1.8}{1+4\lambda}$.

1. Appliquer la descente de gradient gradDescentParam afin de résoudre numériquement le problème $(P)$, on partira de x0 = Snoise. On pourra choisir nbIter = 50 et lam = 4 (et tau adapté en fonction) pour ce premier test. Représenter et comparer le signal $S_{fin}$ ainsi obtenu avec $S$ et $S_{noise}$.

# un premier test avec lam = 4
lam = 4
nbIter = 50
tau = 1.8/(1+4*lam)

Sfin = gradDescentParam(gradE, lam, Snoise, nbIter, tau)

plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.plot(Snoise, 'orange', label = u'$S_{noise}$')
plt.plot(S, 'black', label = '$S$')
plt.legend()

plt.subplot(1,2,2)
plt.plot(Sfin, 'red', label = u'$S_{fin}$')
plt.plot(S, 'black', label = '$S$')
plt.legend()

<matplotlib.legend.Legend at 0x7fe02039cd90>

1. On s'intéresse à l'influence du paramètre $\lambda$. Conserver les mêmes paramètres et faire varier lam in np.arange(0,50,0.2). Pour chacune de ces valeurs de lam, calculer $S_{fin}$ (comme à la question précédente) puis $N_2 (\lambda) := \| S_{fin} - S \|_2$. Représenter $N_2(\lambda)$ en fonction de $\lambda$.

# on essaie avec plusieurs valeurs de lam

lamVal = np.arange(0,50,0.2)
lamSize = lamVal.size
Norme2 = np.zeros(lamSize)
#Efinal = np.zeros(lamSize)
for i,lam in enumerate(lamVal):
    nbIter = 500
    tau = 1.8/(1+4*lam)
    Sfin = gradDescentParam(gradE, lam, Snoise, nbIter, tau)
    Norme2[i] = np.linalg.norm(S-Sfin,2)
    #Efinal[i] = E(Sfin,lam)/E(Snoise,lam)

plt.figure()
plt.plot(lamVal,Norme2, label = u'$|| S_{fin} - S ||_2$')
plt.xlabel('$\lambda$')
plt.legend()

<matplotlib.legend.Legend at 0x7fe011047e50>

1. En déduire le choix de paramètre $\lambda$ qui semble optimal dans le cas présent (on "triche" ici : on connaît la solution $S$ qu'on veut retrouver). Pour la valeur du paramètre $\lambda > 0$ ainsi obtenue, représenter et comparer le signal $S_{fin}$ ainsi obtenu avec $S$ et $S_{noise}$. On pourra utiliser la fonction np.argmin.

Il est normal que la valeur de $\lambda$ "optimale" varie dès qu'on réexécute la première cellule de l'exercice :

noise = np.random.normal(0, 0.1, N) Snoise = S + noise

ind = np.argmin(Norme2)
lamOpt = lamVal[ind]
errOpt = np.min(Norme2)
print('lamOpt = '+str(lamOpt)+u' et dans ce cas ||S - Snoise||_2 = '+str(errOpt))

lam = lamOpt
tau = 1.8/(1+4*lam)

SfinE2 = gradDescentParam(gradE, lam, Snoise, 500, tau)

plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.plot(Snoise, 'orange', label = u'$S_{noise}$')
plt.plot(S, 'black', label = '$S$')
plt.legend()

plt.subplot(1,2,2)
plt.plot(SfinE2, 'red', label = u'$S_{fin}$')
plt.plot(S, 'black', label = '$S$')
plt.legend()

lamOpt = 4.4 et dans ce cas ||S - Snoise||_2 = 0.6725222655921174

<matplotlib.legend.Legend at 0x7fe010fd0ca0>

Exercice 5.- Modèle variationnel pour débruiter une image

On reprend notre image initiale et on moyenne les 3 canaux RVB afin d'obtenir une image en niveaux de gris stockée dans le np.ndarray U de taille n,N, ici $n=N=256$. On génère ensuite une version bruitée en ajoutant à chaque pixel une v.a. suivant une loi normale centrée de variance $\sigma^2$ avec $\sigma = 0.1$. On obtient le np.ndarray Unoise de même taille que U.

n,N,a = img0.shape # n : nb de lignes et N : nb de colonnes

plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.imshow(img0)
plt.title('hibiscus')

U = np.sum(img0, axis=2)/(3*255)
plt.subplot(1,3,2)
plt.imshow(U, cmap='gray')
plt.title('U')

noise = np.random.normal(0, 0.05, size=(n,N)) # sigma = 0.1
Unoise = U + noise

plt.subplot(1,3,3)
plt.imshow(Unoise, cmap='gray')
plt.title('Unoise')

Text(0.5, 1.0, 'Unoise')

Reprenons la démarche suivie pour aboutir au modèle $1D$ (début de l'exercice $1$). On se donne à présent $v : [0,1] \times [0,1] \rightarrow \mathbb{R}$ qu'on suppose être une version bruitée de notre image $u : [0,1] \times [0,1] \rightarrow \mathbb{R}$ qu'on ne connaît pas et qu'on souhaiterait "retrouver". Pour cela, on va à nouveau essayer de différencier (de façon quantitative) une image bruitée d'une image non bruitée en se basant sur les oscillations/variations : pour $q \geq 1$, $$\| \nabla v \|_{q}^q = \int_{[0,1] \times [0,1]} | \nabla v |^q = \int_{[0,1] \times [0,1]} \left(| \partial_x v |^2 + | \partial_y v|^2 \right)^\frac{q}{2} \quad \text{sera plus grand pour } v \text{ que pour } u \: .$$ Du point de vue du signal discret, on peut remplacer $\partial_x v (x_i, y_j)$ par le taux de variation entre deux valeurs successives : supposons qu'on connaît

$$V = \left[ v(x_i, y_j) \right]_{\substack{i = 0 \ldots N-1 \\ j = 0 \ldots n-1}} \in {\rm M}_{N,n}(\mathbb{R}) \quad \text{avec} \quad 0 = x_0 < x_1 < \ldots < x_{N-1} = 1 \quad \text{et} \quad 0 = y_0 < y_1 < \ldots < y_{n-1} = 1$$

on peut remplacer $\partial_x v$ par la matrice des différences selon les abscisses $\delta^x V \in {\rm M}_{n,N-1}(\mathbb{R})$ :

$$\begin{align*} \delta^x V & = \left[ \begin{array}{cccc} v(x_1, y_0) - v(x_0, y_0) & v(x_2, y_0) - v(x_1, y_0) & \ldots & v(x_{N-1}, y_0) - v(x_{N-2}, y_0) \\ v(x_1, y_1) - v(x_0, y_1) & v(x_2, y_1) - v(x_1, y_1) & \ldots & v(x_{N-1}, y_1) - v(x_{N-2}, y_1) \\ \vdots & \vdots & & \vdots \\ v(x_1, y_{n-1}) - v(x_0, y_{n-1}) & v(x_2, y_{n-1}) - v(x_1, y_{n-1}) & \ldots & v(x_{N-1}, y_{n-1}) - v(x_{N-2} y_{n-1}) \\ \end{array} \right] \\ & = \left[ \begin{array}{cccc} V_{0,1} - V_{0,0} & V_{0,2} - V_{0,1} & \ldots & V_{0,{N-1}} - V_{0,{N-2}} \\ V_{1,1} - V_{1,0} & V_{1,2} - V_{1,1} & \ldots & V_{1,{N-1}} - V_{1,{N-2}} \\ \vdots & \vdots & & \vdots \\ V_{n-1,1} - V_{n-1,0} & V_{n-1,2} - V_{n-1,1} & \ldots & V_{n-1,N-1} - V_{n-1,N-2} \\ \end{array} \right] = \left[ \begin{array}{c|c|c|c} V_{\cdot, 1} - V_{\cdot, 0} & V_{\cdot, 2} - V_{\cdot, 1} & \ldots & V_{\cdot,{N-1}} - V_{\cdot,{N-2}} \\ \end{array} \right] \end{align*}$$

On a utilisé les notations $V_{i, \cdot}$ pour désigner la ligne $i$ de $V$ et $V_{\cdot, j}$ pour désigner la colonne $j$ de $V$.

On peut de même remplacer $\partial_y V$ par la matrice des différences selon les ordonnées $\delta^y V \in {\rm M}_{n-1,N}(\mathbb{R})$ :

$$\delta^y V = \left[ \begin{array}{cccc} V_{1,0} - V_{0,0} & V_{1,1} - V_{0,1} & \ldots & V_{1,N-1} - V_{0,N-1} \\ V_{2,0} - V_{1,0} & V_{2,1} - V_{1,1} & \ldots & V_{2,N-1} - V_{1,N-1} \\ \vdots & \vdots & & \vdots \\ V_{n-1,0} - V_{n-2,0} & V_{n-1,1} - V_{n-2,1} & \ldots & V_{{n-1},N-1} - V_{{n-2},N-1} \\ \end{array} \right] = \left[ \begin{array}{c} V_{1,\cdot} - V_{0,\cdot} \\ \hline V_{2,\cdot} - V_{1,\cdot} \\ \hline \vdots \\ \hline V_{n-1,\cdot} - V_{n-2,\cdot} \\ \end{array} \right]$$

On peut alors remplacer $\| \nabla v \|_{q}^q$ par

$$\begin{align*} J(V) & = \sum_{i,j} \left( | (\delta^x V)_{i,j} |^2 + | (\delta^y V)_{i,j}|^2 \right)^\frac{q}{2} \\ & = \sum_{i=1}^{N-1} \sum_{j=1}^{n-1} \left( | V_{i,j} - V_{i-1,j}|^2 + | V_{i,j} - V_{i,j-1}|^2 \right)^\frac{q}{2} + \sum_{j = 1}^{n-1} | V_{0,j} - V_{0,j-1}|^q + \sum_{i=1}^{N-1} | V_{i,0} - V_{i-1,0}|^q \: . \end{align*}$$

Pour $q=2$, on peut en déduire facilement/fastidieusement le gradient de $J$, similairement à ce qui a été fait dans le cas d'un signal $1D$.

1. Effectuer le calcul de $\nabla J(V)$ et vérifier qu'on peut alors réécrire, $$\frac{1}{2} \nabla J(V)= - \delta^x \left( \left[  \begin{array}{c|c|c} 0 & & 0 \\ \vdots & (\delta^x V) & \vdots \\ 0 & & 0 \end{array} \right] \right) - \delta^y \left( \left[  \begin{array}{ccc} 0 & \ldots & 0 \\ \hline & (\delta^y V) & \\ \hline 0 & \ldots & 0 \end{array} \right] \right)$$

Soit $k \in \{1, \ldots, N-2 \}$ et $l \in \{1, \ldots, n-2\}$, le terme $V_{k,l}$ apparaît $4$ fois : $2$ fois pour $(i,j) = (k,l)$, pour $(i,j) = (k+1,l)$ et pour $(i,j) = (k,l)$, $(i,j) = (k,l+1)$. On obtient $$\begin{align*} \frac{1}{2}\partial_{k,l} J(V) & = 2 V_{k,l} - V_{k-1,l} - V_{k+1,l} + 2 V_{k,l} - V_{k,l-1} - V_{k,l+1} = 4 V_{k,l} - V_{k-1,l} - V_{k+1,l} - V_{k,l-1} - V_{k,l+1} \\ \frac{1}{2}\partial_{0,l} J(V) & = 3 V_{0,l} - V_{1,l} - V_{0,l-1} - V_{0,l+1} \quad \text{et} \quad \frac{1}{2}\partial_{N-1,l} J(V) = 3 V_{N-1,l} - V_{N-2,l} - V_{N-1,l-1} - V_{N-1,l+1}\\ \frac{1}{2}\partial_{k,0} J(V) & = 3 V_{k,0} - V_{k-1,0} - V_{k+1,0} - V_{k,1} \quad \text{et} \quad \frac{1}{2}\partial_{k,n-1} J(V) = 3 V_{k,n-1} - V_{k-1,n-1} - V_{k+1,n-1} - V_{k,n-2}\\ \frac{1}{2}\partial_{0,0} J(V) & = 2 V_{0,0} - V_{1,0} - V_{0,1}, \quad \frac{1}{2}\partial_{N-1,0} J(V) = 2 V_{N-1,0} - V_{N-2,0} - V_{N-1,1}, \quad \frac{1}{2}\partial_{0,n-1} J(V) = 2 V_{0,0} - V_{1,n-1} - V_{0,n-2} \quad \text{et} \quad \frac{1}{2}\partial_{N-1,n-1} J(V) = 2 V_{N-1,n-1} - V_{N-2,n-1} - V_{N-1,n-2}. \\ \end{align*}$$

De même que dans le cas $1D$, on peut alors réécrire,

$$\begin{align*} \nabla J(V) & = \left[  \begin{array}{c|c|c} 0 & & 0 \\ \vdots & - \delta^x (\delta^x V) & \vdots \\ 0 & & 0 \end{array} \right] + \left[  \begin{array}{ccc} 0 & \ldots & 0 \\ \hline & - \delta^y (\delta^y V) & \\ \hline 0 & \ldots & 0 \end{array} \right] + \left[  \begin{array}{c|c|c} & & \\ - (\delta^x V)_{\cdot, 0} & 0_{{\rm M}_{n,N-2}} & (\delta^x V)_{\cdot, N-1} \\ & & \end{array} \right] + \left[  \begin{array}{ccc} & - (\delta^y V)_{0,\cdot} & \\ \hline & 0_{{\rm M}_{n-2,N}} & \\ \hline & (\delta^y V)_{n-1,\cdot} & \end{array} \right] \\ & = - \delta^x \left[  \begin{array}{c|c|c} 0 & & 0 \\ \vdots & (\delta^x V) & \vdots \\ 0 & & 0 \end{array} \right] - \delta^y \left[  \begin{array}{ccc} 0 & \ldots & 0 \\ \hline & (\delta^y V) & \\ \hline 0 & \ldots & 0 \end{array} \right] \end{align*}$$

où par définition $\delta^x (\delta^x V) \in {\rm M}_{n,N-2}(\mathbb{R})$ et $\delta^y (\delta^y V) \in {\rm M}_{n-2,N}(\mathbb{R})$, et $M_{\cdot,j}$ désigne la colonne $j$ de $M$, de même $M_{i,\cdot}$ désigne la ligne $i$ de $M$ d'une matrice $M$.

Cela permet d'implémenter assez facilement le gradient de $J$ et donc celui de $E$.

1. Définir des fonctions deltax et deltay qui prennent en argument un np.ndarray V de taille n,N et renvoient respectivement le np.ndarray $\delta^x V$ de taille n,N-1 et le np.ndarray $\delta^y V$ de taille n-1,N. Grâce à l'écriture ci-dessus du gradient (question $1$), définir une fonction gradEMat(V,lam) qui renvoie le gradient de $$E(V) = \frac{1}{2} \| V - U_{noise} \|_2^2 + \frac{\lambda}{2} J(V) \: .$$

On a $\displaystyle \nabla E (V) = V - U_{noise} + \frac{\lambda}{2} \nabla J(V)$ et on peut utiliser la question $1$ pour l'implémentation.

def deltax(V):
    return V[: , 1:] - V[: , :-1]
    
def deltay(V):
    return V[1: , :] - V[:-1, :]

def gradEMat(V,lam):
    n,N = V.shape
    W1 = np.zeros((n,N+1))
    W2 = np.zeros((n+1,N))
    W1[:, 1:-1] = deltax(V)
    W2[1:-1, :] = deltay(V)
    W = -deltax(W1) - deltay(W2)
    return V - Unoise + lam*W

1. Utiliser la fonction gradDescentParam implémentée précédemment afin d'obtenir une solution approchée Ufin au problème de minimsation de $E$ en partant de x0 = Unoise. On pourra commencer par tester avec $\lambda = 2$, $\displaystyle \tau = \frac{1.8}{1+8\lambda}$ et $100$ itérations.

lam = 2
tau = 1.8/(1+8*lam)
UfinE2 = gradDescentParam(gradEMat, lam, Unoise, 100, tau)

plt.figure(figsize=(15,15))
plt.subplot(3,3,1)
plt.imshow(Unoise, cmap='gray')
plt.subplot(3,3,2)
plt.imshow(UfinE2, cmap='gray')
plt.subplot(3,3,3)
plt.imshow(U, cmap='gray')
plt.subplot(3,3,4)
plt.imshow(Unoise[128:,:128], cmap='gray')
plt.subplot(3,3,5)
plt.imshow(UfinE2[128:,:128], cmap='gray')
plt.subplot(3,3,6)
plt.imshow(U[128:,:128], cmap='gray')
plt.subplot(3,3,7)
plt.imshow(Unoise[100:150,100:150], cmap='gray')
plt.subplot(3,3,8)
plt.imshow(UfinE2[100:150,100:150], cmap='gray')
plt.subplot(3,3,9)
plt.imshow(U[100:150,100:150], cmap='gray')

<matplotlib.image.AxesImage at 0x7fe011000730>

1. Représenter l'image Ufin obtenue à la question précédente et la comparer avec U et Unoise, on pourra comparer quelques fenêtres de zoom.

# voir question précédente

1. Jouer sur les différents paramètres et comparer les résultats obtenus.

# À vous de jouer

Exercice 6.- L'effet du régularisateur ${\rm L}^2$

On va essayer d'étudier l'effet du choix de "régularisation" $J = J_2$. On remarque tout d'abord que la méthode de débruitage mise en oeuvre dans l'exercice précédent produit un floutage des contours présents dans l'image. On va le vérifier dans les questions $1$ et $2$ suivantes.

1. Appliquer la descente de gradient à l'énergie modifée $$\tilde{E} (V) = \frac{1}{2} \| V - U \|_2^2 + \frac{\lambda}{2} J(V)$$ où on a remplacé Unoise par U dans le terme d'attache aux données afin d'observer uniquement l'effet de $J_2$. On partira de x0 = U directement, en conservant les paramètres $\lambda = 2$, $\displaystyle \tau = \frac{1.8}{1+8\lambda}$. Comparer les résultats après $10$ itérations et après $100$ itérations avec l'image U.

On modifie la fonction gradDescentParampour passer le terme d'attache aux données en paramètre.

# descente de gradient avec parametre lam et xdata pour le terme d'attache aux données
def gradDescentParam2(gradfun, lam, xdata, x0, nbIter, tau):
    xk = x0
    for k in range(nbIter):
        xk = xk - tau*gradfun(xk, lam, xdata)
    return xk

def gradEMat2(V,lam, xdata):
    n,N = V.shape
    W1 = np.zeros((n,N+1))
    W2 = np.zeros((n+1,N))
    W1[:, 1:-1] = deltax(V)
    W2[1:-1, :] = deltay(V)
    W = -deltax(W1) - deltay(W2)
    return V - xdata + lam*W

lam = 2
tau = 1.8/(1+8*lam)
Ufin2 = gradDescentParam2(gradEMat2, lam,U, U, 10, tau)
Ufin3 = gradDescentParam2(gradEMat2, lam,U, U, 100, tau)
plt.figure(figsize=(15,15))
plt.subplot(3,3,1)
plt.imshow(Ufin2, cmap='gray')
plt.title('Après 10 itérations')
plt.subplot(3,3,2)
plt.imshow(Ufin3, cmap='gray')
plt.title('Après 100 itérations')
plt.subplot(3,3,3)
plt.imshow(U, cmap='gray')
plt.title('Image initiale U')
plt.subplot(3,3,4)
plt.imshow(Ufin2[128:,:128], cmap='gray')
plt.subplot(3,3,5)
plt.imshow(Ufin3[128:,:128], cmap='gray')
plt.subplot(3,3,6)
plt.imshow(U[128:,:128], cmap='gray')
plt.subplot(3,3,7)
plt.imshow(Ufin2[100:150,100:150], cmap='gray')
plt.subplot(3,3,8)
plt.imshow(Ufin3[100:150,100:150], cmap='gray')
plt.subplot(3,3,9)
plt.imshow(U[100:150,100:150], cmap='gray')

<matplotlib.image.AxesImage at 0x7fe00af748b0>

On observe que même en l'absence de bruit dans le term d'attache aux données et en partant de l'image elle-même, la minimisation de $E = E_2$ introduit un floutage des contours présents dans l'image.

1. Effectuer la même expérience qu'à la question $1$, en remplaçant U par l'image plus simple Inb définie ci-dessous.

Inb = np.zeros((256,256))
for i in range(100,150):
    for j in range(50,200):
        Inb[i,j] = 1

lam = 2
tau = 1.8/(1+8*lam)
Inbfin2 = gradDescentParam2(gradEMat2, lam, Inb, Inb, 10, tau)
Inbfin3 = gradDescentParam2(gradEMat2, lam, Inb, Inb, 100, tau)
plt.figure(figsize=(15,10))
plt.subplot(2,3,1)
plt.imshow(Inbfin2, cmap='gray')
plt.title('Après 10 itérations')
plt.subplot(2,3,2)
plt.imshow(Inbfin3, cmap='gray')
plt.title('Après 100 itérations')
plt.subplot(2,3,3)
plt.imshow(Inb, cmap='gray')
plt.title('Image initiale Inb')
plt.subplot(2,3,4)
plt.imshow(Inbfin2[75:125,25:75], cmap='gray')
plt.subplot(2,3,5)
plt.imshow(Inbfin3[75:125,25:75], cmap='gray')
plt.subplot(2,3,6)
plt.imshow(Inb[75:125,25:75], cmap='gray')

<matplotlib.image.AxesImage at 0x7fe00ad1ef10>

1. On revient dans cette question à un signal 1D très simple : reprendre à nouveau la question 1 en remplaçant l'image Upar le signal Sgeom défini ci-dessous. Représenter le signal Sgeom et les signaux S10 et S50 obtenus après $10$ itérations et $50$ itérations de descente de gradient.

Plus précisément, on reprendra le gradient de E avec paramètre $\lambda$ : gradE(V, lam) défini à la question $1$ de l'exercice $4$ et on remplacera Snoise par Sgeom dans le terme d'attache aux données. On effectuera ensuite la descente de gradient grâce à gradDescentParam (exercice $4$ question $2$) en prenant x0 = Sgeom, $\lambda = 10$, $\tau = \displaystyle \frac{1.8}{1+4\lambda}$ et nbIt = 10 puis nbIt = 50.

t01 = np.linspace(0,1,256)
Sgeom = (t01>0.5)**np.ones(256)
plt.figure()
plt.plot(t01,Sgeom, 'black', label = u'$S_{geom}$')

# gradient de E 1D avec parametre lam et terme d'attache Sgeom
def gradEgeom(V,lam):
    N = V.size
    W = np.zeros(N)
    W[1:-1] = - delta(delta(V))
    W[0] = V[0] - V[1]
    W[N-1] = V[N-1] - V[N-2]
    return V - Sgeom + lam*W

lam = 10
tau = 1.8/(1+4*lam)

S10 = gradDescentParam(gradEgeom, lam, Sgeom, 10, tau)
S50 = gradDescentParam(gradEgeom, lam, Sgeom, 50, tau)

plt.plot(t01, S10, 'blue', label = u'$S_{10}$')
plt.plot(t01, S50, 'red', label = u'$S_{50}$')
plt.legend()

<matplotlib.legend.Legend at 0x7fe00ac1feb0>

On observe que la transition $0$-$1$ est régularisée lorsqu'on minimise $E$ pour $q = 2$. rappelons à présent que dans l'exercice $1$, on avait proposé un terme de régularisation de la forme $$J_q(S) = \| \delta S \|_q^q$$ et on avait observé que plus $q$ était proche de $1$, plus $J_q(S) - J_q(S_{noise})$ était grand ce qui donne envie de considérer $q = 1$. Cependant, le choix $q = 2$ était également naturel car simple à mettre en oeuvre : le gradient de $J_2$ est assez simple à exprimer, ce qui permet d'utiliser une descente de gradient comme on l'a fait jusqu'ici.

1. Calculer et comparer $J_q$ pour $q = 1$ et $q = 2$ pour chacun des trois signaux Sgeom, S10 et S50 représentés à la question $3$. En ce qui concerne les valeurs obtenues pour $J_1$, justifiez les résultats obtenus par un calcul théorique.

def J(V, q):
    return np.linalg.norm(delta(V),q)**q

for q in [1,2]:
    print('J'+str(q)+'(Sgeom) = '+str(J(Sgeom,q)))
    print('J'+str(q)+'(S10) = '+str(J(S10,q)))
    print('J'+str(q)+'(S50) = '+str(J(S50,q)))

J1(Sgeom) = 1.0
J1(S10) = 1.0
J1(S50) = 1.0
J2(Sgeom) = 1.0
J2(S10) = 0.1114420763720354
J2(S50) = 0.08123556609447082

En ce qui concerne $J_1$, dès que le signal est croissant au sens où pour tout $k \in \{1, \ldots, N-1\}$, $V_k \geq V_{k-1}$, on a $$J_1(V) = \sum_{k=1}^{N-1} | V_k - V_{k-1} | = \sum_{k=1}^{N-1} ( V_k - V_{k-1} ) = V_{N-1} - V_0 \: ,$$ indépendemment du fait que la transition entre la valeur initiale et la valeur finale soit régulière ou brutale. C'est très différent en ce qui concerne $J_2$ qui prend des valeurs bien plus petites lorsque la transition se fait plus régulièrement.

On observe ainsi que contrairement à $J_2$, le terme $J_1$ ne tend pas à régulariser les transitions brutales, et on pourrait vérifier que de même $J_1$ ne tend pas à flouter les contours présents dans une image (ces contours correspondent aux lignes de transition brutale entre deux teintes de gris différentes).

Exercice 7.- Modèle ${\rm TV}$ - ${\rm L}^2$ régularisé

On revient au choix $q = 1$. Notons tout d'abord une obstruction importante : $J_1$ n'est pas différentiable et la mise en oeuvre de la descente de gradient n'est plus possible. Pour $V \in \mathbb{R}^N$, $$J_1 (V) = \sum_{i=1}^{N-1} | V_i - V_{i-1} | \:$$ n'est pas différentiable en tout point $V$ tel que $\exists k$, $V_k = V_{k-1}$. Afin de contourner cette obstruction, on propose de remplacer $|t|$ par $\sqrt{\epsilon^2 + t^2}$, et ainsi considérer pour régulariser le signal $S_{noise}$ : $$E_\epsilon (V) = \frac{1}{2} \| V - S_{noise} \|_2^2 + \lambda J_\epsilon (V) \quad \text{avec} \quad J_\epsilon(V) = \sum_{i=1}^{N-1} \sqrt{ \epsilon^2 + | V_i - V_{i-1} |^2 } = \sum_{k=i}^{N-2} \sqrt{ \epsilon^2 + | (\delta V)_{i} |^2 } \: ,$$ ou encore pour régulariser l'image $U_{noise}$ : $$E_\epsilon (V) = \frac{1}{2} \| V - U_{noise} \|_2^2 + \lambda J_\epsilon (V) \quad \text{avec} \quad J_\epsilon(V) = \sum_{i,j} \sqrt{\epsilon^2 + | (\delta^x V)_{i,j} |^2 + | (\delta^y V)_{i,j}|^2 } \: .$$

1. On définit pour $\epsilon > 0$, $h_\epsilon : \mathbb{R} \rightarrow \mathbb{R}_+$, $t \mapsto \sqrt{\epsilon^2 + t^2}$. Représenter le graphe de $h_\epsilon$ sur $[-2,2]$ pour $\epsilon \in \{1, 0.1, 0.01\}$. Calculer $h_\epsilon^\prime$ et représenter également son graphe.

def h(t,eps):
    return np.sqrt(eps**2 + t**2)

def hprime(t,eps):
    return t/h(t,eps)

plt.figure(figsize=(20,5))
plt.subplot(1,2,1)
t = np.linspace(-2,2,41)
for eps in [1, 0.1, 0.01]:
    plt.plot(t, h(t,eps), label = u'$\epsilon = $'+str(eps))
plt.plot(t,np.abs(t),'+', label = u'$t \mapsto |t|$')
plt.legend()
plt.title('$h_\epsilon$')

plt.subplot(1,2,2)
t = np.linspace(-2,2,41)
for eps in [1, 0.1, 0.01]:
    plt.plot(t, hprime(t,eps), label = u'$\epsilon = $'+str(eps))
plt.legend()
plt.title('$h_\epsilon^\prime$')

Text(0.5, 1.0, '$h_\\epsilon^\\prime$')

1. Soit $i,k \in \{0, 1, \ldots, N-1\}$, vérifier que $$\partial_k \left[ h_\epsilon(V_i - V_{i-1}) \right] = \left\lbrace \begin{array}{ccc} h_\epsilon^\prime (V_k - V_{k-1}) & \text{si} & i = k \\ -h_\epsilon^\prime (V_{k+1} - V_{k}) & \text{si} & i = k+1 \\ 0 & \text{sinon} & \end{array}\right. \: .$$ En déduire que pour $0< k < N-1$, $$\partial_k J_\epsilon (V) = h_\epsilon^\prime (V_k - V_{k-1}) - h_\epsilon^\prime (V_{k+1} - V_{k}) \quad \text{et} \quad \partial_0 J_\epsilon(V) = -h_\epsilon^\prime (V_1 - V_0) \quad \text{et} \quad \partial_{N-1} J_\epsilon(V) = h_\epsilon^\prime ( V_{N-1} - V_{N-2}) \: .$$

1. Fixer $\epsilon = 0.01$ et implémenter une fonction gradEeps(V) qui prend en argument un np.array V et retourne $\nabla E_\epsilon (V)$ où $E_\epsilon(V) = \frac{1}{2}\| V - S_{noise} \|_2^2 + \lambda J_\epsilon(V)$.

On pourra utiliser la question $2$ pour réécrire $$\nabla J_\epsilon(V) = -\delta W \quad \text{où} \quad W = \left[ 0 \right. ,\underbrace{h_\epsilon^\prime ( (\delta V)_1 ) , h_\epsilon^\prime ( (\delta V)_2 ) , \ldots , h_\epsilon^\prime ( (\delta V)_{N-2} ) }_{\in \mathbb{R}^{N-2}}, \left. 0 \right] \: ,$$ où $\delta V$ a été défini à l'exercice $1$.

eps = 0.01

def gradJeps(V):
    N = V.size
    W = np.zeros(N+1)
    W[1:-1] = hprime(delta(V), eps)
    return - delta(W)


def gradEeps(V,lam):
    return V - Snoise + lam*gradJeps(V)

1. Appliquer la descente de gradient avec $\lambda = 0.2$, $\tau = \displaystyle \frac{1.8}{1 + 4 \frac{\lambda}{\epsilon}}$, en partant de x0 = Snoise et en effectuant $200$ itérations. Comparer avec le signal obtenu en minimisant $E = E_2$ dans l'exercice $4$.

lam = 0.2
tau = 1.8/(1+4*lam/eps)

Sfin = gradDescentParam(gradEeps, lam, Snoise, 200, tau)

plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.plot(Snoise, 'orange', label = u'$S_{noise}$')
plt.plot(S, 'black', label = '$S$')
plt.legend()

plt.subplot(1,2,2)
plt.plot(Sfin, 'red', label = u'$S_{fin}$')
plt.plot(S, 'black', label = '$S$')
plt.legend()

<matplotlib.legend.Legend at 0x7fe00aa5e7c0>

Afin de comparer au débruitage obtenu avec l'énergie $E_2$, on effectue l'étude sur le paramètre $\lambda$ "optimal" analogue à la question $6$ de l'exercice $4$.

# on essaie avec plusieurs valeurs de lam

print(r'On rappelle que dans le cas du débruitage E2 on avait')
print('lamOpt = '+str(lamOpt)+u' et dans ce cas ||S - Snoise||_2 = '+str(errOpt))


lamVal = np.arange(0,2,0.01)
lamSize = lamVal.size
Norme2 = np.zeros(lamSize)
for i,lam in enumerate(lamVal):
    nbIter = 200
    tau = 1.8/(1+4*lam/eps)
    Sfin = gradDescentParam(gradEeps, lam, Snoise, nbIter, tau)
    Norme2[i] = np.linalg.norm(S-Sfin,2)

plt.figure()
plt.plot(lamVal,Norme2, label = u'$|| S_{fin} - S ||_2$')
plt.xlabel('$\lambda$')
plt.legend()

ind = np.argmin(Norme2)
lamOpt = lamVal[ind]
errOpt = np.min(Norme2)
print(r'Dans le cas du débruitage Eeps on obtient')
print('lamOpt = '+str(lamOpt)+u' et dans ce cas ||S - Snoise||_2 = '+str(errOpt))

lam = lamOpt
tau = 1.8/(1+4*lam/eps)

Sfin = gradDescentParam(gradEeps, lam, Snoise, 200, tau)

plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.plot(Snoise, 'orange', label = u'$S_{noise}$')
plt.plot(S, 'black', label = '$S$')
plt.legend()

plt.subplot(1,3,2)
plt.plot(Sfin, 'red', label = u'$S_{fin}$')
plt.plot(S, 'black', label = '$S$')
plt.title('débruitage $E_\epsilon$')
plt.legend()

plt.subplot(1,3,3)
plt.plot(SfinE2, 'red', label = u'$S_{fin}$')
plt.plot(S, 'black', label = '$S$')
plt.legend()
plt.title('débruitage $E_2$')

On rappelle que dans le cas du débruitage E2 on avait
lamOpt = 4.4 et dans ce cas ||S - Snoise||_2 = 0.6725222655921174
Dans le cas du débruitage Eeps on obtient
lamOpt = 0.2 et dans ce cas ||S - Snoise||_2 = 0.6419488352060425

Text(0.5, 1.0, 'débruitage $E_2$')

1. On va terminer ce travail en menant à bien la minimisation de $\displaystyle E_\epsilon(V) = \frac{1}{2} \| V - U_{noise} \|_2^2 + \lambda J_\epsilon(V)$ avec $\displaystyle h_\epsilon(t_1,t_2) = \sqrt{\epsilon^2 + t_1^2 + t_2^2}$ et $$\begin{align*} J_\epsilon (V) & = \sum_{i,j} h_\epsilon \left((\delta^x V)_{i,j} , (\delta^y V)_{i,j} \right) \\ & = \sum_{i=1}^{N-1} \sum_{j=1}^{n-1} h_\epsilon \left(  V_{i,j} - V_{i,j-1} , V_{i,j} - V_{i-1,j} \right) + \sum_{j = 1}^{n-1} h_\epsilon \left(  V_{0,j} - V_{0,j-1}, 0 \right) + \sum_{i=1}^{N-1} h_\epsilon \left( 0, V_{i,0} - V_{i-1,0} \right) \: . \end{align*}$$ On conserve $\epsilon = 0.01$ et on admettra que le gradient de $E_\epsilon$ peut-être implémenté comme suit par la fonction gradETV(V,lam). Appliquer la descente de gradient afin de minimiser $E_\epsilon$, on prendra $\lambda = 0.05$, $\tau = \displaystyle \frac{1.8}{1+8\frac{\lambda}{\epsilon}}$, X0 = Unoise et $300$ itérations. Comparer avec les résultats obtenus à la fin de l'exercice $5$. Vérifier qu'on obtient des contours plus nets dans l'image.

eps = 0.01

def gradETV(V, lam):
    n,N = V.shape
    W1 = np.zeros((n,N))
    W2 = np.zeros((n,N))
    dx = deltax(V); dy = deltay(V)
    W1[:, :-1] = dx
    W2[:-1, :] = dy
    Neps = np.sqrt(eps**2 + W1**2 + W2**2)
    #
    W3 = np.zeros((n,N+1))
    W4 = np.zeros((n+1,N))
    W3[:, 1:-1] = dx/Neps[:, :-1]
    W4[1:-1, :] = dy/Neps[:-1, :]
    W = -deltax(W3) - deltay(W4)
    return V - Unoise + lam*W


lam = 0.05
tau = 1.8/(1+8*lam/eps)

UfinTV = gradDescentParam(gradETV, lam, Unoise, 300, tau)
plt.figure(figsize=(15,15))
plt.subplot(3,3,1)
plt.imshow(Unoise, cmap='gray')
plt.subplot(3,3,2)
plt.imshow(UfinTV, cmap='gray')
plt.subplot(3,3,3)
plt.imshow(U, cmap='gray')
plt.subplot(3,3,4)
plt.imshow(Unoise[128:,:128], cmap='gray')
plt.subplot(3,3,5)
plt.imshow(UfinTV[128:,:128], cmap='gray')
plt.subplot(3,3,6)
plt.imshow(U[128:,:128], cmap='gray')
plt.subplot(3,3,7)
plt.imshow(Unoise[100:150,100:150], cmap='gray')
plt.subplot(3,3,8)
plt.imshow(UfinTV[100:150,100:150], cmap='gray')
plt.subplot(3,3,9)
plt.imshow(U[100:150,100:150], cmap='gray')

<matplotlib.image.AxesImage at 0x7fe0103f8a60>

1. Représenter (par une image) $\nabla E_\epsilon(U_{noise})$ et $\nabla E_2(U_{noise})$ en utilisant les paramètres ayant servis pour les calculs numériques dans la question précédente/dans la question $3$ de l'exercice $5$.

plt.figure(figsize=(16,8))
plt.subplot(1,2,1)
plt.imshow(gradETV(Unoise, 0.05), cmap ='gray')
plt.title(r'$\nabla E_\epsilon(U_{noise})$')
plt.subplot(1,2,2)
plt.imshow(gradEMat(Unoise, 2), cmap ='gray')
plt.title(r'$\nabla E_2(U_{noise})$')

plt.figure(figsize=(16,8))
plt.subplot(1,2,1)
plt.imshow(U-UfinTV, cmap ='gray')
plt.title(r'$U- U_{fin}$ pour $E_\epsilon$')
plt.subplot(1,2,2)
plt.imshow(U-UfinE2, cmap ='gray')
plt.title(r'$U- U_{fin}$ pour $E_2$')

Text(0.5, 1.0, '$U- U_{fin}$ pour $E_2$')