Bildrekonstruktion mittels Compressive Sensing

Einführung

In diesem Lab wird die Rekonstruktion eines dünn besetzten Bildes aus einer Reihe paralleler Projektionen unter Verwendung von Compressive Sensing demonstriert. Compressive Sensing ist eine Technik zur effizienten Erfassung und Rekonstruktion von Signalen, die in einem bestimmten Bereich dünn besetzt sind. In diesem Fall interessieren wir uns für die Rekonstruktion eines 2D-Bildes aus einer geringen Anzahl von Projektionen, die entlang unterschiedlicher Winkel aufgenommen wurden. Wir werden die Leistung von L1- und L2-Penalisierungsmethoden für diese Aufgabe vergleichen.

Tipps für die VM

Nachdem der VM-Start abgeschlossen ist, klicken Sie in der oberen linken Ecke, um zur Registerkarte Notebook zu wechseln und Jupyter Notebook für die Übung zu öffnen.

Manchmal müssen Sie einige Sekunden warten, bis Jupyter Notebook vollständig geladen ist. Die Validierung von Vorgängen kann aufgrund von Einschränkungen in Jupyter Notebook nicht automatisiert werden.

Wenn Sie bei der Lernphase Probleme haben, können Sie Labby gerne fragen. Geben Sie nach der Sitzung Feedback, und wir werden das Problem für Sie prompt beheben.

Skills Graph

Bibliotheken importieren

In diesem Schritt werden wir die erforderlichen Bibliotheken importieren.

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

Hilfsfunktionen definieren

In diesem Schritt werden wir Hilfsfunktionen definieren, um die Tomographie-Designmatrix zu generieren und synthetische binäre Daten zu erstellen.

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))

Daten generieren

In diesem Schritt werden wir synthetische binäre Daten und Projektionen generieren.

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)

Bild mit L2-Penalisierung rekonstruieren

In diesem Schritt werden wir das Bild unter Verwendung von L2- (Ridge-)Penalisierung rekonstruieren.

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

Bild mit L1-Penalisierung rekonstruieren

In diesem Schritt werden wir das Bild unter Verwendung von L1- (Lasso-)Penalisierung rekonstruieren.

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

Ergebnisse visualisieren

In diesem Schritt werden wir das ursprüngliche Bild und die rekonstruierten Bilder unter Verwendung von L2- und L1-Penalisierung visualisieren.

plt.figure(figsize=(8, 3.3))
plt.subplot(131)
plt.imshow(data, cmap=plt.cm.gray, interpolation="nearest")
plt.axis("off")
plt.title("Original Image")
plt.subplot(132)
plt.imshow(rec_l2, cmap=plt.cm.gray, interpolation="nearest")
plt.title("L2 Penalization")
plt.axis("off")
plt.subplot(133)
plt.imshow(rec_l1, cmap=plt.cm.gray, interpolation="nearest")
plt.title("L1 Penalization")
plt.axis("off")

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

plt.show()

Zusammenfassung

In diesem Lab haben wir die Verwendung von Compressive Sensing zur Rekonstruktion eines 2D-Bilds aus einer geringen Anzahl von Projekionen, die entlang unterschiedlicher Winkel aufgenommen wurden, demonstriert. Wir haben die Leistung von L1- und L2-Penalisierungsmethoden verglichen und festgestellt, dass die L1-Penalisierung ein besseres Ergebnis mit weniger Artefakten liefert.