当前位置: 首页>>技术教程>>正文


sklearn例程:压缩感知:使用L1先验(Lasso)进行层析成像重建

压缩感知

什么是压缩感知呢?

压缩感知

示例说明:使用L1先验(Lasso)进行层析成像重建

此示例介绍了从一组沿不同角度获取的平行投影来重建图像的方法。这种数据集一般是通过CT检查获得的。CT即计算机断层扫描或计算层析成像

在样本上没有任何先验信息的情况下,重建图像所需的投影数约为图像线性大小l(width x height)的量级(以像素为单位)。为简单起见,我们在这里考虑一个稀疏图像,其中只有图像中对象边界上的像素具有非零值。不过,大多数图像在不同的基础上都是稀疏的,例如Haar小波。因为数据的稀疏性,只能得到l/7投影,因此像做图像重建有必要使用样本上可用的先验信息(稀疏度):这就是压缩感知

层析成像投影操作是线性变换。除了对应于线性回归的数据保真项之外,我们还对图像的L1范数进行了惩罚,以考虑其稀疏性。由此产生的优化问题称为Lasso,它对应的scikit-learn类是sklearn.linear_model.Lasso,使用坐标下降算法求解实现。重要的是,与这里使用的投影运算符相比,该实现在稀疏矩阵上的计算效率更高。

即使将噪声添加到投影中,使用L1惩罚进行的重建也会产生零误差的结果(所有像素均成功标记为0或1)。相比之下,L2惩罚(sklearn.linear_model.Ridge)会为像素产生大量错误标记。也就是说,L2与L1惩罚相比,重建的图像上会观察到明显的伪影。

代码实现[Python]


# -*- coding: utf-8 -*- 

print(__doc__)

# Author: Emmanuelle Gouillart 
# License: BSD 3 clause

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.
    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.) ** 2 + (y - l / 2.) ** 2  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)

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

# 使用L1惩罚(Lasso)重构
# 最佳alpha值,采用交叉验证(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分9.761秒。
运行代码输出的图片内容如下:

Compressive sensing: tomography reconstruction with L1 prior (Lasso)

源码下载

参考资料

本文由《纯净天空》出品。文章地址: https://vimsky.com/article/4481.html,未经允许,请勿转载。