压缩感知:使用 L1 先验(Lasso)进行断层扫描重建#

此示例展示了如何从一组沿不同角度获取的平行投影重建图像。这种数据集是在**计算机断层扫描**(CT)中获取的。

在没有任何关于样本的先验信息的情况下,重建图像所需的投影数量与图像的线性大小 l(以像素为单位)的数量级相同。为简单起见,我们在这里考虑一个稀疏图像,其中只有物体边界上的像素具有非零值。例如,此类数据可能对应于细胞材料。但是请注意,大多数图像在不同的基(例如 Haar 小波)中是稀疏的。仅获取了 l/7 个投影,因此有必要使用样本上可用的先验信息(其稀疏性):这是一个**压缩感知**的示例。

断层扫描投影操作是一种线性变换。除了对应于线性回归的数据保真度项之外,我们还惩罚图像的 L1 范数以考虑其稀疏性。由此产生的优化问题称为 Lasso。我们使用 Lasso 类,它使用坐标下降算法。重要的是,与此处使用的投影算子相比,此实现对于稀疏矩阵的计算效率更高。

即使投影中添加了噪声,使用 L1 惩罚进行的重建也能给出零误差的结果(所有像素都成功标记为 0 或 1)。相比之下,L2 惩罚(Ridge)会产生大量像素标记错误。与 L1 惩罚相反,在重建图像上观察到重要的伪影。特别要注意的是,圆形伪影将角落中的像素分开,这些像素对投影的贡献少于中心圆盘。

original image, L2 penalization, L1 penalization
# Author: Emmanuelle Gouillart <[email protected]>
# License: BSD 3 clause

import matplotlib.pyplot as plt
import numpy as np
from scipy import ndimage, sparse

from sklearn.linear_model import Lasso, Ridge


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):
    """Compute the tomography design matrix.

    Parameters
    ----------

    l_x : int
        linear size of image array

    n_dir : int
        number of angles at which projections are acquired.

    Returns
    -------
    p : sparse matrix of shape (n_dir l_x, l_x**2)
    """
    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():
    """Synthetic binary 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))


# Generate synthetic images, and 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 with L2 (Ridge) penalization
rgr_ridge = Ridge(alpha=0.2)
rgr_ridge.fit(proj_operator, proj.ravel())
rec_l2 = rgr_ridge.coef_.reshape(l, l)

# Reconstruction with L1 (Lasso) penalization
# the best value of alpha was determined using cross validation
# with LassoCV
rgr_lasso = Lasso(alpha=0.001)
rgr_lasso.fit(proj_operator, proj.ravel())
rec_l1 = rgr_lasso.coef_.reshape(l, l)

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

**脚本总运行时间:**(0 分 10.145 秒)

相关示例

硬币图像上的结构化 Ward 层次聚类的演示

硬币图像上的结构化 Ward 层次聚类的演示

用于图像分割的谱聚类

用于图像分割的谱聚类

基于 L1 的稀疏信号模型

基于 L1 的稀疏信号模型

层次聚类:结构化与非结构化 Ward

层次聚类:结构化与非结构化 Ward

由 Sphinx-Gallery 生成的图库