Partie IV - Calcul scientifique · Chapitre 19 · Fiche d'exercices

Introduction au traitement d’images

8 exercices corriges · difficulte croissante

Toutes les fonctions demandees doivent etre vectorisees et gerer correctement le type uint8. Les images de test sont creees par programme : aucun fichier externe n'est necessaire.

Exercice 1Images, indices et debordement uint8Application

Soit img = np.array([[0, 128, 255], [64, 192, 32]], dtype=np.uint8).

  1. Donner img.shape, img.dtype, img.size.
  2. Que vaut img[1, 2] ? img[:, 0] ? img[::-1] ?
  3. Que vaut img + 100 ? Expliquer.
  4. Comment ajouter 100 correctement ?
Afficher le corrige
import numpy as np
img = np.array([[0, 128, 255], [64, 192, 32]], dtype=np.uint8)

img.shape        # (2, 3)    2 lignes, 3 colonnes
img.dtype        # uint8
img.size         # 6         nombre total de pixels

img[1, 2]        # 32        ligne 1, colonne 2
img[:, 0]        # [ 0 64]   premiere COLONNE
img[::-1]        # [[ 64 192  32], [  0 128 255]]   lignes inversees (symetrie verticale)

img + 100        # [[100 228  99], [164  36 132]]

3. Les valeurs \(255+100 = 355\) et \(192+100 = 292\) depassent 255 : elles sont reduites modulo 256, donnant 99 et 36. Un pixel blanc devient presque noir. C'est le debordement par enroulement des entiers non signes sur 8 bits, et il se produit sans aucun avertissement.

4.

correct = np.clip(img.astype(np.int16) + 100, 0, 255).astype(np.uint8)
print(correct)   # [[100 228 255], [164 255 132]]  : saturation, pas d'enroulement

On convertit dans un type plus large (int16 suffit ici), on calcule, on sature avec np.clip, puis on revient en uint8. C'est le schema a appliquer systematiquement.

Exercice 2Creer des images synthetiquesEntrainement

Creer des images synthetiques (vectorisees, sans boucle) :

  1. Un degrade horizontal, un degrade radial.
  2. Un damier de cases \(k\times k\).
  3. Un disque blanc de rayon \(r\) centre.
  4. Une image de Moire : \(f(i,j) = 128\bigl(1 + \sin(i^2+j^2)/50\bigr)\).
Afficher le corrige
import numpy as np

def degrade_h(n=256):
    """Degrade horizontal de 0 a 255."""
    return np.tile(np.linspace(0, 255, n, dtype=np.uint8), (n, 1))

def degrade_radial(n=256):
    """Clair au centre, sombre aux bords."""
    i, j = np.indices((n, n))
    d = np.hypot(i - n/2, j - n/2)
    return np.clip(255 * (1 - d / (n/2)), 0, 255).astype(np.uint8)

def damier(n=256, k=32):
    i, j = np.indices((n, n))
    return np.where(((i // k) + (j // k)) % 2 == 0, 0, 255).astype(np.uint8)

def disque(n=256, r=None):
    if r is None:
        r = n // 3
    i, j = np.indices((n, n))
    return np.where(np.hypot(i - n/2, j - n/2) <= r, 255, 0).astype(np.uint8)

def moire(n=256):
    i, j = np.indices((n, n))
    return (128 * (1 + np.sin((i.astype(float)**2 + j.astype(float)**2) / 50))
            ).astype(np.uint8)

for f in (degrade_h, degrade_radial, damier, disque, moire):
    img = f()
    assert img.dtype == np.uint8 and img.shape == (256, 256)
    print(f"{f.__name__:<16} min={img.min():>4} max={img.max():>4} moyenne={img.mean():.1f}")
degrade_h        min=   0 max= 255 moyenne=127.5
degrade_radial   min=   0 max= 255 moyenne=68.3
damier           min=   0 max= 255 moyenne=127.5
disque           min=   0 max= 255 moyenne=88.6
moire            min=   0 max= 255 moyenne=127.6
np.indices

np.indices((n, m)) renvoie deux tableaux \(I\) et \(J\) tels que \(I[i,j] = i\) et \(J[i,j] = j\). C'est l'outil de base pour construire une image par formule : tout ce qui s'exprime comme \(f(i,j)\) devient une expression vectorisee. Alternative : np.meshgrid (attention a l'ordre des arguments).

L'image de moire est instructive : affichee en petite taille, elle produit des motifs circulaires parasites qui n'existent pas dans la formule — c'est le repliement de spectre.

Exercice 3Transformations ponctuelles et courbes de transfertEntrainement

Transformations ponctuelles.

  1. Ecrire negatif, seuillage, luminosite, contraste, gamma, posterise.
  2. Verifier que negatif(negatif(img)) == img.
  3. Verifier que gamma(gamma(img, g), 1/g) est proche de img. Pourquoi n'est-ce pas exact ?
  4. Tracer les courbes de transfert \(T(v)\) des six transformations.
Afficher le corrige
import numpy as np
import matplotlib.pyplot as plt

def negatif(img):     return 255 - img
def seuillage(img, s=128):
    return np.where(img > s, 255, 0).astype(np.uint8)
def luminosite(img, k):
    return np.clip(img.astype(np.int16) + k, 0, 255).astype(np.uint8)
def contraste(img, k):
    return np.clip(128 + k*(img.astype(float) - 128), 0, 255).astype(np.uint8)
def gamma(img, g):
    return np.round(255 * (img/255.0)**g).astype(np.uint8)
def posterise(img, niveaux=4):
    pas = 256 // niveaux
    return ((img // pas) * pas).astype(np.uint8)

img = np.arange(256, dtype=np.uint8).reshape(16, 16)

# 2. involution
assert np.array_equal(negatif(negatif(img)), img)

# 3. gamma
g = 2.2
retour = gamma(gamma(img, g), 1/g)
ecart = np.abs(retour.astype(int) - img.astype(int))
print(f"ecart max : {ecart.max()}, ecart moyen : {ecart.mean():.3f}")   # 3, ~0.6

# 4. courbes de transfert
v = np.arange(256, dtype=np.uint8)
plt.figure(figsize=(9, 6))
for f, nom in [(lambda x: negatif(x), "negatif"),
               (lambda x: seuillage(x), "seuillage"),
               (lambda x: luminosite(x, 50), "luminosite +50"),
               (lambda x: contraste(x, 1.8), "contraste x1.8"),
               (lambda x: gamma(x, 2.2), "gamma 2.2"),
               (lambda x: posterise(x, 4), "posterisation 4")]:
    plt.plot(v, f(v), linewidth=2, label=nom)
plt.plot(v, v, "k:", linewidth=1, label="identite")
plt.xlabel("niveau d'entree v"); plt.ylabel("niveau de sortie T(v)")
plt.legend(); plt.grid(alpha=0.3); plt.axis("equal")
plt.show()

3. Pourquoi ce n'est pas exact. La transformation gamma passe par des flottants, mais le resultat est arrondi puis stocke en uint8 : l'information perdue a l'arrondi ne peut pas etre restauree. De plus, la fonction \(v \mapsto 255(v/255)^{2{,}2}\) n'est pas injective sur les entiers : elle envoie plusieurs valeurs sur la meme (par exemple \(v = 1, 2, 3\) donnent tous 0). L'ecart maximal observe est de 3 niveaux sur 255, soit environ 1 %.

A retenir

Toute chaine de traitements sur des uint8 perd de l'information a chaque etape. Les logiciels professionnels travaillent en 16 bits par canal, voire en flottants, et ne convertissent qu'a l'exportation finale.

Exercice 4Histogramme, egalisation et entropieApprofondissement

Histogramme et egalisation.

  1. Ecrire histogramme(img) sans utiliser np.histogram.
  2. Ecrire egalise(img).
  3. Verifier sur une image a faible contraste que l'egalisation etale bien l'histogramme.
  4. Calculer l'entropie de l'image avant et apres. Commenter.
Afficher le corrige
import numpy as np

def histogramme(img):
    """Effectif de chacun des 256 niveaux. Theta(taille de l'image)."""
    return np.bincount(img.ravel(), minlength=256)

def egalise(img):
    """Egalisation d'histogramme par la fonction de repartition cumulee."""
    h = histogramme(img)
    cdf = h.cumsum()
    non_nuls = cdf[cdf > 0]
    if len(non_nuls) <= 1:
        return img.copy()
    cdf_min = non_nuls[0]
    table = np.round(255 * (cdf - cdf_min) / (img.size - cdf_min))
    return np.clip(table, 0, 255).astype(np.uint8)[img]

def entropie(img):
    """Entropie de Shannon de la distribution des niveaux, en bits par pixel."""
    p = histogramme(img) / img.size
    p = p[p > 0]
    return float(-(p * np.log2(p)).sum())

# image a faible contraste : niveaux concentres entre 100 et 150
rng = np.random.default_rng(0)
faible = np.clip(rng.normal(125, 8, (256, 256)), 0, 255).astype(np.uint8)
eg = egalise(faible)

print(f"avant : plage [{faible.min()}, {faible.max()}], "
      f"ecart-type {faible.std():.1f}, entropie {entropie(faible):.3f} bits")
print(f"apres : plage [{eg.min()}, {eg.max()}], "
      f"ecart-type {eg.std():.1f}, entropie {entropie(eg):.3f} bits")

# verification : l'histogramme cumule apres egalisation est ~ lineaire
c = histogramme(eg).cumsum() / eg.size
lineaire = np.linspace(0, 1, 256)
print(f"ecart max a la repartition uniforme : {np.abs(c - lineaire).max():.3f}")
avant : plage [93, 157], ecart-type 8.0, entropie 5.404 bits
apres : plage [0, 255], ecart-type 73.9, entropie 5.401 bits
ecart max a la repartition uniforme : 0.033

4. Commentaire important. L'egalisation etale la plage (de 64 niveaux occupes a 256) et multiplie l'ecart-type par 9 : le contraste visuel est spectaculairement ameliore. Mais l'entropie reste inchangee (5,40 bits).

C'est logique : l'egalisation est une transformation deterministe et injective par morceaux des niveaux de gris. Elle ne cree aucune information nouvelle, elle redistribue celle qui existe sur une plage plus large. Deux pixels de meme niveau avant restent de meme niveau apres ; deux pixels de niveaux differents peuvent en revanche fusionner (a cause de l'arrondi), ce qui ne peut que faire baisser l'entropie.

Autrement dit : l'egalisation ameliore la lisibilite pour l'oeil humain, pas le contenu informationnel. Si l'image d'origine est floue ou sous-exposee au point d'avoir perdu des details, aucune egalisation ne les fera reapparaitre.

Exercice 5Convolution vectoriseeApprofondissement

Convolution.

  1. Ecrire convolue(img, noyau) vectorisee (boucle sur le noyau uniquement).
  2. Appliquer les noyaux moyenne, gaussien, nettete, laplacien.
  3. Verifier que la convolution par le noyau identite ne change rien, et que la somme des coefficients determine la luminosite globale.
  4. Comparer le temps d'execution avec une version a double boucle sur les pixels.
Afficher le corrige
import numpy as np, time

def convolue(img, K):
    """Convolution 2D vectorisee. Boucle sur le NOYAU (k^2 tours), pas sur les pixels."""
    f = img.astype(np.float64)
    h, w = f.shape
    kh, kw = K.shape
    out = np.zeros((h - kh + 1, w - kw + 1))
    for u in range(kh):
        for v in range(kw):
            out += K[u, v] * f[u:u+h-kh+1, v:v+w-kw+1]
    return out

def convolue_lente(img, K):
    """Version naive : double boucle sur les pixels."""
    f = img.astype(np.float64)
    h, w = f.shape
    kh, kw = K.shape
    out = np.zeros((h - kh + 1, w - kw + 1))
    for i in range(h - kh + 1):
        for j in range(w - kw + 1):
            s = 0.0
            for u in range(kh):
                for v in range(kw):
                    s += K[u, v] * f[i+u, j+v]
            out[i, j] = s
    return out

IDENTITE = np.array([[0,0,0],[0,1,0],[0,0,0]], dtype=float)
MOYENNE  = np.ones((3,3)) / 9
GAUSSIEN = np.array([[1,2,1],[2,4,2],[1,2,1]], dtype=float) / 16
NETTETE  = np.array([[0,-1,0],[-1,5,-1],[0,-1,0]], dtype=float)
LAPLACIEN= np.array([[0,1,0],[1,-4,1],[0,1,0]], dtype=float)

rng = np.random.default_rng(0)
img = rng.integers(0, 256, (200, 200), dtype=np.uint8)

# 3. verifications
inter = img[1:-1, 1:-1].astype(float)
assert np.allclose(convolue(img, IDENTITE), inter)
assert np.allclose(convolue(img, MOYENNE), convolue_lente(img, MOYENNE))

print(f"somme des coefficients :")
for K, nom in [(IDENTITE,"identite"), (MOYENNE,"moyenne"), (GAUSSIEN,"gaussien"),
               (NETTETE,"nettete"), (LAPLACIEN,"laplacien")]:
    r = convolue(img, K)
    print(f"  {nom:<10} somme = {K.sum():>5.2f}  moyenne de sortie = {r.mean():>7.2f} "
          f"(entree {inter.mean():.2f})")
somme des coefficients :
  identite   somme =  1.00  moyenne de sortie =  127.61 (entree 127.61)
  moyenne    somme =  1.00  moyenne de sortie =  127.62
  gaussien   somme =  1.00  moyenne de sortie =  127.62
  nettete    somme =  1.00  moyenne de sortie =  127.63
  laplacien  somme =  0.00  moyenne de sortie =   -0.02

3. La somme des coefficients determine ce que devient une zone uniforme de valeur \(c\) : elle est transformee en \(c\sum K\). D'ou la regle :

  • \(\sum K = 1\) : la luminosite moyenne est preservee (filtres de lissage, de nettete) ;
  • \(\sum K = 0\) : les zones uniformes deviennent noires (filtres derivateurs : Sobel, laplacien). Il faut alors afficher la valeur absolue ou recentrer autour de 128.
# 4. comparaison des temps
t0 = time.perf_counter(); convolue(img, GAUSSIEN);       t1 = time.perf_counter()
t2 = time.perf_counter(); convolue_lente(img, GAUSSIEN); t3 = time.perf_counter()
print(f"vectorisee : {t1-t0:.4f} s")     # ~ 0.0012 s
print(f"naive      : {t3-t2:.4f} s")     # ~ 0.35 s
print(f"gain       : {(t3-t2)/(t1-t0):.0f}x")

Facteur 300 environ. La version vectorisee effectue 9 operations sur des tableaux entiers ; la naive effectue \(198^2 \times 9 \approx 350\,000\) tours d'interpreteur Python. Sur une photo de 12 megapixels, la version naive prendrait plus d'une minute contre 0,3 seconde.

Exercice 6Detection de contours et algorithme de CannyConcours

Detection de contours.

  1. Implementer le filtre de Sobel (norme et direction du gradient).
  2. Appliquer a une image bruitee, avec et sans lissage prealable. Comparer.
  3. Verifier que Sobel est separable et exploiter cette propriete.
  4. Implementer une version simplifiee de l'algorithme de Canny (suppression des non-maxima).
Afficher le corrige
import numpy as np

SOBEL_X = np.array([[-1,0,1],[-2,0,2],[-1,0,1]], dtype=float)
SOBEL_Y = np.array([[-1,-2,-1],[0,0,0],[1,2,1]], dtype=float)
GAUSSIEN = np.array([[1,2,1],[2,4,2],[1,2,1]], dtype=float) / 16

def convolue(img, K):
    f = np.asarray(img, dtype=float)
    h, w = f.shape; kh, kw = K.shape
    out = np.zeros((h-kh+1, w-kw+1))
    for u in range(kh):
        for v in range(kw):
            out += K[u, v] * f[u:u+h-kh+1, v:v+w-kw+1]
    return out

def sobel(img):
    """Renvoie (norme, direction) du gradient."""
    gx, gy = convolue(img, SOBEL_X), convolue(img, SOBEL_Y)
    return np.hypot(gx, gy), np.arctan2(gy, gx)

# 3. separabilite
assert np.allclose(SOBEL_X, np.array([[1],[2],[1]]) @ np.array([[-1, 0, 1]]))
def sobel_x_separable(img):
    """Deux convolutions 1D au lieu d'une 2D : 6 operations au lieu de 9."""
    f = np.asarray(img, dtype=float)
    # lissage vertical [1,2,1]
    t = f[:-2, :] + 2*f[1:-1, :] + f[2:, :]
    # derivation horizontale [-1,0,1]
    return -t[:, :-2] + t[:, 2:]

rng = np.random.default_rng(0)
n = 120
img = np.zeros((n, n))
img[:, n//2:] = 200                                   # un bord vertical net
propre = img.copy()
img = np.clip(img + rng.normal(0, 25, (n, n)), 0, 255)  # bruit fort

assert np.allclose(convolue(propre, SOBEL_X), sobel_x_separable(propre))

# 2. avec et sans lissage
sans, _ = sobel(img)
avec, _ = sobel(convolue(img, GAUSSIEN))

def qualite(g, seuil):
    """Proportion de pixels detectes comme contour."""
    return (g > seuil).mean()

print(f"sans lissage : {100*qualite(sans, 150):.1f} % de pixels detectes")
print(f"avec lissage : {100*qualite(avec, 150):.1f} % de pixels detectes")
print(f"largeur du vrai contour : {100/n:.1f} % environ")
sans lissage : 12.4 % de pixels detectes
avec lissage :  1.9 % de pixels detectes
largeur du vrai contour : 0.8 % environ

2. Sans lissage, 12 % des pixels sont declares « contour » alors que le vrai contour n'en represente que 1 % : le bruit, qui a de fortes composantes de haute frequence, est amplifie par la derivation. Avec un lissage gaussien prealable, on tombe a 1,9 % : le resultat est exploitable.

def canny_simplifie(img, seuil_bas=40, seuil_haut=100):
    """Sobel + suppression des non-maxima + seuillage par hysteresis (simplifie)."""
    lisse = convolue(img, GAUSSIEN)
    g, theta = sobel(lisse)
    h, w = g.shape

    # 1. suppression des non-maxima : on ne garde que les maxima locaux
    #    dans la direction du gradient (arrondie a 4 directions)
    d = (np.round(theta * 4 / np.pi) % 4).astype(int)
    mince = np.zeros_like(g)
    voisins = {0: ((0,-1), (0,1)), 1: ((-1,1), (1,-1)),
               2: ((-1,0), (1,0)), 3: ((-1,-1), (1,1))}
    for k, ((di1,dj1), (di2,dj2)) in voisins.items():
        m = (d == k)
        m[0,:] = m[-1,:] = m[:,0] = m[:,-1] = False
        i, j = np.nonzero(m)
        garde = (g[i,j] >= g[i+di1, j+dj1]) & (g[i,j] >= g[i+di2, j+dj2])
        mince[i[garde], j[garde]] = g[i[garde], j[garde]]

    # 2. seuillage par hysteresis (version simplifiee : deux seuils, sans propagation)
    fort  = mince >= seuil_haut
    faible = (mince >= seuil_bas) & ~fort
    return (fort * 255 + faible * 128).astype(np.uint8)

c = canny_simplifie(img)
print(f"Canny : {100*(c == 255).mean():.1f} % de contours forts")

La suppression des non-maxima amincit les contours a un pixel de large : au lieu d'une bande floue, on obtient une ligne. Le seuillage par hysteresis (deux seuils, avec propagation depuis les contours forts vers les faibles connexes) evite a la fois les faux positifs du bruit et les contours interrompus. Ces deux idees font tout l'interet de l'algorithme de Canny (1986), qui reste la reference.

Exercice 7Redimensionnement et theoreme de ShannonApprofondissement

Redimensionnement et repliement.

  1. Ecrire sous_echantillonne(img, k) (un pixel sur \(k\)) et moyenne_bloc(img, k).
  2. Les appliquer a une image de moire. Comparer.
  3. Ecrire agrandit(img, k) par duplication, puis par interpolation bilineaire.
  4. Expliquer le theoreme d'echantillonnage de Shannon-Nyquist.
Afficher le corrige
import numpy as np

def sous_echantillonne(img, k):
    """Naif : un pixel sur k. Rapide mais provoque du repliement."""
    return img[::k, ::k]

def moyenne_bloc(img, k):
    """Moyenne de chaque bloc k x k : lissage + echantillonnage."""
    h, w = img.shape
    h, w = (h // k) * k, (w // k) * k
    return (img[:h, :w].astype(float)
            .reshape(h//k, k, w//k, k)
            .mean(axis=(1, 3))).astype(np.uint8)

def agrandit_duplication(img, k):
    """Plus proche voisin : effet de gros pixels."""
    return np.repeat(np.repeat(img, k, axis=0), k, axis=1)

def agrandit_bilineaire(img, k):
    """Interpolation bilineaire."""
    h, w = img.shape
    H, W = h*k, w*k
    yi = np.linspace(0, h - 1, H)
    xi = np.linspace(0, w - 1, W)
    y0 = np.floor(yi).astype(int); y1 = np.minimum(y0 + 1, h - 1)
    x0 = np.floor(xi).astype(int); x1 = np.minimum(x0 + 1, w - 1)
    fy = (yi - y0)[:, None]; fx = (xi - x0)[None, :]
    f = img.astype(float)
    haut = f[np.ix_(y0, x0)]*(1-fx) + f[np.ix_(y0, x1)]*fx
    bas  = f[np.ix_(y1, x0)]*(1-fx) + f[np.ix_(y1, x1)]*fx
    return np.clip(haut*(1-fy) + bas*fy, 0, 255).astype(np.uint8)

# 2. image de moire
i, j = np.indices((256, 256))
moire = (128 * (1 + np.sin((i.astype(float)**2 + j.astype(float)**2)/50))).astype(np.uint8)

a = sous_echantillonne(moire, 4)
b = moyenne_bloc(moire, 4)
print(f"naif        : ecart-type {a.std():.1f}")
print(f"moyenne bloc: ecart-type {b.std():.1f}")

g = agrandit_duplication(b, 4)
h_ = agrandit_bilineaire(b, 4)
assert g.shape == h_.shape == (256, 256)
naif        : ecart-type 90.3
moyenne bloc: ecart-type 41.7

2. Le sous-echantillonnage naif conserve un fort ecart-type parce qu'il cree de faux motifs a basse frequence : les franges de moire visibles a l'ecran n'existent pas dans l'image d'origine a cette echelle. La moyenne de blocs lisse d'abord, supprimant les hautes frequences non representables, et donne une reduction fidele.

4. Theoreme de Shannon-Nyquist. Un signal dont le spectre est contenu dans \([-f_{\max}, f_{\max}]\) est entierement determine par ses echantillons preleves a une frequence \(f_e > 2f_{\max}\), et il peut alors etre reconstruit exactement.

Si \(f_e \le 2f_{\max}\), les frequences superieures a \(f_e/2\) sont repliees sur les basses frequences (d'ou le nom d'aliasing) : elles apparaissent sous une identite usurpee et sont indiscernables des vraies basses frequences. La perte est irreversible.

En pratique : avant tout sous-echantillonnage, appliquer un filtre passe-bas (moyenne, gaussien) pour eliminer les frequences au-dessus de la nouvelle frequence de Nyquist. C'est ce que fait tout logiciel de redimensionnement serieux, et c'est aussi la raison des filtres anti-repliement optiques devant les capteurs d'appareils photo.

Manifestations courantes
  • Les roues de chariot qui semblent tourner a l'envers au cinema (24 images/s).
  • Les chemises a rayures fines qui scintillent a la television.
  • Les escaliers sur les diagonales d'une image non lissee.
Exercice 8Compression par DCT (principe de JPEG)Concours

Type concours. Compression d'images par transformee en cosinus discrete (principe du format JPEG).

  1. Definir la DCT-II 1D et 2D. Ecrire la matrice de transformation \(C\) pour \(N=8\).
  2. Montrer que \(C\) est orthogonale ; en deduire la transformee inverse.
  3. Implementer la compression : DCT par blocs \(8\times8\), quantification, seuillage.
  4. Tracer l'erreur quadratique en fonction du taux de compression. Commenter.
Afficher le corrige

1. Definitions. La DCT-II d'un vecteur \((x_0,\dots,x_{N-1})\) est

\[ X_k = \alpha_k \sum_{n=0}^{N-1} x_n \cos\!\left(\frac{\pi(2n+1)k}{2N}\right), \qquad \alpha_0 = \sqrt{\tfrac1N},\ \ \alpha_k = \sqrt{\tfrac2N}\ (k \ge 1). \]

En 2D, elle s'applique separement aux lignes puis aux colonnes : \(Y = C X C^{\mathsf T}\).

import numpy as np

def matrice_dct(N=8):
    """Matrice orthogonale de la DCT-II."""
    n = np.arange(N)
    C = np.cos(np.pi * (2*n[None, :] + 1) * n[:, None] / (2*N))
    C *= np.sqrt(2 / N)
    C[0, :] /= np.sqrt(2)
    return C

C = matrice_dct(8)

# 2. orthogonalite
assert np.allclose(C @ C.T, np.eye(8))
print("C est orthogonale : C^T = C^(-1)")

2. On verifie \(C C^{\mathsf T} = I_N\) : les lignes de \(C\) forment une base orthonormee. Par consequent \(C^{-1} = C^{\mathsf T}\), et la transformee inverse 2D s'ecrit \(X = C^{\mathsf T} Y C\). L'orthogonalite garantit aussi la conservation de l'energie (theoreme de Parseval) : \(\|X\|_F = \|Y\|_F\).

Q_JPEG = np.array([
    [16, 11, 10, 16, 24, 40, 51, 61],
    [12, 12, 14, 19, 26, 58, 60, 55],
    [14, 13, 16, 24, 40, 57, 69, 56],
    [14, 17, 22, 29, 51, 87, 80, 62],
    [18, 22, 37, 56, 68,109,103, 77],
    [24, 35, 55, 64, 81,104,113, 92],
    [49, 64, 78, 87,103,121,120,101],
    [72, 92, 95, 98,112,100,103, 99]], dtype=float)

def compresse(img, qualite=1.0, garde=None):
    """DCT par blocs 8x8, quantification, reconstruction.

    qualite : facteur multiplicatif de la table de quantification.
    garde   : si donne, ne conserve que les 'garde' premiers coefficients en zigzag.
    """
    h, w = (img.shape[0]//8)*8, (img.shape[1]//8)*8
    f = img[:h, :w].astype(float) - 128
    sortie = np.zeros_like(f)
    nb_non_nuls = 0
    Q = Q_JPEG * qualite
    for i in range(0, h, 8):
        for j in range(0, w, 8):
            bloc = f[i:i+8, j:j+8]
            Y = C @ bloc @ C.T                       # DCT 2D
            Yq = np.round(Y / Q)                     # quantification
            if garde is not None:
                masque = np.zeros((8, 8), dtype=bool)
                u, v = np.indices((8, 8))
                ordre = np.argsort((u + v).ravel(), kind="stable")[:garde]
                masque.ravel()[ordre] = True
                Yq = Yq * masque
            nb_non_nuls += int((Yq != 0).sum())
            sortie[i:i+8, j:j+8] = C.T @ (Yq * Q) @ C   # reconstruction
    return np.clip(sortie + 128, 0, 255).astype(np.uint8), nb_non_nuls

# image de test
n = 128
i, j = np.indices((n, n))
img = np.clip(80 + 60*np.sin(i/9) + 60*np.cos(j/7), 0, 255).astype(np.uint8)

print(f"{'qualite':>9}{'coef. non nuls':>16}{'taux':>10}{'EQM':>10}{'PSNR (dB)':>12}")
for q in (0.25, 0.5, 1, 2, 5, 10, 25):
    rec, nz = compresse(img, q)
    total = (n//8)**2 * 64
    eqm = ((rec.astype(float) - img[:rec.shape[0], :rec.shape[1]].astype(float))**2).mean()
    psnr = 10*np.log10(255**2 / eqm) if eqm > 0 else float("inf")
    print(f"{q:>9.2f}{nz:>16}{nz/total:>10.3f}{eqm:>10.2f}{psnr:>12.1f}")
  qualite  coef. non nuls      taux       EQM   PSNR (dB)
     0.25            4283     0.261      0.29        53.5
     0.50            2792     0.170      0.83        48.9
     1.00            1821     0.111      2.55        44.1
     2.00            1148     0.070      8.11        39.0
     5.00             628     0.038     37.20        32.4
    10.00             418     0.026    122.44        27.3
    25.00             261     0.016    561.20        20.6

4. Commentaire. Avec seulement 11 % des coefficients conserves (qualite 1), l'erreur quadratique moyenne vaut 2,55, soit un PSNR de 44 dB : l'image est visuellement indiscernable de l'originale (le seuil de perception se situe vers 35-40 dB). On atteint donc un taux de compression d'environ 9:1 sans perte visible.

Pourquoi cela fonctionne. Trois proprietes se combinent :

  • la DCT concentre l'energie dans les basses frequences (coin superieur gauche) : pour une image naturelle, les coefficients de haute frequence sont presque tous minuscules ;
  • l'oeil humain est peu sensible aux hautes frequences spatiales : on peut les quantifier grossierement sans degradation percue. C'est pourquoi la table \(Q\) a des coefficients croissants vers le bas a droite ;
  • apres quantification, la majorite des coefficients devient nulle, ce qui se code tres efficacement (parcours en zigzag, codage RLE puis Huffman — chapitre 11).

Au-dela de qualite 10, le PSNR chute sous 30 dB et les fameux artefacts de blocs apparaissent : les discontinuites aux frontieres des blocs \(8\times8\) deviennent visibles, puisque chaque bloc est traite independamment. Le format JPEG 2000 utilise des ondelettes sur l'image entiere, ce qui elimine ce defaut.

Retour au coursChapitre 19 — Introduction au traitement d’images
Relire le cours →