Введение
В этом лабораторном задании демонстрируется реконструкция разреженного изображения из набора параллельных проекций с использованием сжатие чувствительности. Сжатие чувствительности - это метод эффективного получения и реконструкции сигналов, которые являются разреженными в некоторой области. В данном случае мы заинтересованы в реконструкции 2D-изображения из небольшого числа проекций, полученных по разным углам. Мы сравним производительность методов штрафования L1 и L2 для этой задачи.
Советы по работе с ВМ
После запуска ВМ нажмите в левом верхнем углу, чтобы переключиться на вкладку Notebook и получить доступ к Jupyter Notebook для практики.
Иногда вам может потребоваться подождать несколько секунд, пока Jupyter Notebook не загрузится полностью. Валидация операций не может быть автоматизирована из-за ограничений Jupyter Notebook.
Если вы сталкиваетесь с проблемами во время обучения, не стесняйтесь обращаться к Labby. Оставьте отзыв после занятия, и мы оперативно решим проблему для вас.
Импорт библиотек
В этом шаге мы импортируем необходимые библиотеки.
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
Определение вспомогательных функций
В этом шаге мы определим вспомогательные функции для генерации матрицы проекций томографии и создания синтетических бинарных данных.
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))
Генерация данных
В этом шаге мы сгенерируем синтетические бинарные данные и проекции.
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)
реконструкция изображения с использованием штрафования L2
В этом шаге мы будем реконструировать изображение с использованием штрафования L2 (Ridge).
rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)
реконструкция изображения с использованием штрафования L1
В этом шаге мы будем реконструировать изображение с использованием штрафования L1 (Lasso).
rgr_lasso = Lasso(alpha=0.001)
rgr_lasso.fit(proj_operator, proj.ravel())
rec_l1 = rgr_lasso.coef_.reshape(l, l)
Визуализация результатов
В этом шаге мы визуализируем исходное изображение и реконструированные изображения с использованием штрафования L2 и L1.
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()
Резюме
В этом лабораторном занятии мы продемонстрировали применение сжатие данных для реконструкции 2D-изображения из небольшого числа проекций, полученных при разных углах. Мы сравнили производительность методов штрафования L1 и L2 и наблюдали, что штрафование L1 дает лучший результат с меньшим количеством артефактов.