Reconstruction d'image par compression de données

Introduction

Ce laboratoire démontre la reconstruction d'une image creuse à partir d'un ensemble de projections parallèles en utilisant la compression par capteurs. La compression par capteurs est une technique pour acquérir et reconstruire efficacement des signaux qui sont creux dans un certain domaine. Dans ce cas, nous sommes intéressés à reconstruire une image 2D à partir d'un nombre réduit de projections acquises selon différents angles. Nous comparerons les performances des méthodes de pénalisation L1 et L2 pour cette tâche.

Conseils sur la machine virtuelle

Une fois le démarrage de la machine virtuelle terminé, cliquez dans le coin supérieur gauche pour basculer vers l'onglet Carnet de notes pour accéder au carnet Jupyter pour la pratique.

Parfois, vous devrez peut-être attendre quelques secondes pour que le carnet Jupyter ait fini de charger. La validation des opérations ne peut pas être automatisée en raison des limitations du carnet Jupyter.

Si vous rencontrez des problèmes pendant l'apprentissage, n'hésitez pas à demander à Labby. Donnez votre feedback après la session, et nous résoudrons rapidement le problème pour vous.

Skills Graph

Importation des bibliothèques

Dans cette étape, nous allons importer les bibliothèques nécessaires.

import numpy as np
from scipy import sparse
from scipy import ndimage
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
import matplotlib.pyplot as plt

Définition de fonctions d'aide

Dans cette étape, nous allons définir des fonctions d'aide pour générer la matrice de conception de la tomographie et créer des données binaires synthétiques.

def _weights(x, dx=1, orig=0):
    x = np.ravel(x)
    floor_x = np.floor((x - orig) / dx).astype(np.int64)
    alpha = (x - orig - floor_x * dx) / dx
    return np.hstack((floor_x, floor_x + 1)), np.hstack((1 - alpha, alpha))

def _generate_center_coordinates(l_x):
    X, Y = np.mgrid[:l_x, :l_x].astype(np.float64)
    center = l_x / 2.0
    X += 0.5 - center
    Y += 0.5 - center
    return X, Y

def build_projection_operator(l_x, n_dir):
    X, Y = _generate_center_coordinates(l_x)
    angles = np.linspace(0, np.pi, n_dir, endpoint=False)
    data_inds, weights, camera_inds = [], [], []
    data_unravel_indices = np.arange(l_x**2)
    data_unravel_indices = np.hstack((data_unravel_indices, data_unravel_indices))
    for i, angle in enumerate(angles):
        Xrot = np.cos(angle) * X - np.sin(angle) * Y
        inds, w = _weights(Xrot, dx=1, orig=X.min())
        mask = np.logical_and(inds >= 0, inds < l_x)
        weights += list(w[mask])
        camera_inds += list(inds[mask] + i * l_x)
        data_inds += list(data_unravel_indices[mask])
    proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))
    return proj_operator

def generate_synthetic_data():
    rs = np.random.RandomState(0)
    n_pts = 36
    x, y = np.ogrid[0:l, 0:l]
    mask_outer = (x - l / 2.0) ** 2 + (y - l / 2.0) ** 2 < (l / 2.0) ** 2
    mask = np.zeros((l, l))
    points = l * rs.rand(2, n_pts)
    mask[(points[0]).astype(int), (points[1]).astype(int)] = 1
    mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)
    res = np.logical_and(mask > mask.mean(), mask_outer)
    return np.logical_xor(res, ndimage.binary_erosion(res))

Générer des données

Dans cette étape, nous allons générer des données binaires synthétiques et des projections.

l = 128
proj_operator = build_projection_operator(l, l // 7)
data = generate_synthetic_data()
proj = proj_operator @ data.ravel()[:, np.newaxis]
proj += 0.15 * np.random.randn(*proj.shape)

Reconstruction d'une image avec pénalisation L2

Dans cette étape, nous allons reconstruire l'image en utilisant la pénalisation L2 (Ridge).

rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)

Reconstruction d'une image avec pénalisation L1

Dans cette étape, nous allons reconstruire l'image en utilisant la pénalisation L1 (Lasso).

rgr_lasso = Lasso(alpha=0.001)
rgr_lasso.fit(proj_operator, proj.ravel())
rec_l1 = rgr_lasso.coef_.reshape(l, l)

Visualiser les résultats

Dans cette étape, nous allons visualiser l'image d'origine et les images reconstruites en utilisant la pénalisation L2 et L1.

plt.figure(figsize=(8, 3.3))
plt.subplot(131)
plt.imshow(data, cmap=plt.cm.gray, interpolation="nearest")
plt.axis("off")
plt.title("Image d'origine")
plt.subplot(132)
plt.imshow(rec_l2, cmap=plt.cm.gray, interpolation="nearest")
plt.title("Pénalisation L2")
plt.axis("off")
plt.subplot(133)
plt.imshow(rec_l1, cmap=plt.cm.gray, interpolation="nearest")
plt.title("Pénalisation L1")
plt.axis("off")

plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)

plt.show()

Sommaire

Dans ce laboratoire, nous avons démontré l'utilisation de la compression de données pour reconstruire une image 2D à partir d'un nombre réduit de projections acquises selon différents angles. Nous avons comparé les performances des méthodes de pénalisation L1 et L2 et avons observé que la pénalisation L1 produit un résultat meilleur avec moins d'artefacts.