.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/applications/plot_tomography_l1_reconstruction.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_applications_plot_tomography_l1_reconstruction.py: ====================================================================== 压缩感知：具有L1先验（Lasso）的断层扫描重建 ====================================================================== 本示例展示了如何从一组沿不同角度获取的平行投影中重建图像。这类数据集是在**计算机断层扫描**（CT）中获取的。 在没有任何关于样本的先验信息的情况下，重建图像所需的投影数量与图像的线性尺寸 ``l`` （以像素为单位）成正比。为简单起见，我们在此考虑一个稀疏图像，其中只有物体边界上的像素具有非零值。这类数据可以例如对应于一种细胞材料。然而请注意，大多数图像在不同的基底中是稀疏的，例如Haar小波。仅获取了 ``l/7`` 的投影，因此有必要使用样本的可用先验信息（其稀疏性）：这是**压缩感知**的一个例子。 断层扫描投影操作是一个线性变换。除了对应于线性回归的数据保真项外，我们还对图像的L1范数进行惩罚以考虑其稀疏性。由此产生的优化问题称为 :ref:`lasso` 。我们使用类:class:`~sklearn.linear_model.Lasso` ，它使用坐标下降算法。重要的是，这种实现比这里使用的投影算子在稀疏矩阵上更具计算效率。 具有L1惩罚的重建结果具有零误差（所有像素都成功标记为0或1），即使在投影中添加了噪声。相比之下，L2惩罚（:class:`~sklearn.linear_model.Ridge` ）会导致大量像素标记错误。在重建图像上观察到重要的伪影，这与L1惩罚相反。特别注意将角落像素分开的圆形伪影，这些像素比中央圆盘贡献了更少的投影。 .. GENERATED FROM PYTHON SOURCE LINES 15-122 .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_tomography_l1_reconstruction_001.png :alt: original image, L2 penalization, L1 penalization :srcset: /auto_examples/applications/images/sphx_glr_plot_tomography_l1_reconstruction_001.png :class: sphx-glr-single-img .. code-block:: Python # 作者：scikit-learn 开发者 # SPDX 许可证标识符：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): """计算断层扫描设计矩阵。 Parameters ---------- l_x : int 图像数组的线性尺寸 n_dir : int 获取投影的角度数量 Returns ------- p : 稀疏矩阵，形状为 (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(): """合成二进制数据""" 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（岭）惩罚的重建 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() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 27.757 seconds) .. _sphx_glr_download_auto_examples_applications_plot_tomography_l1_reconstruction.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/scikit-learn/scikit-learn/main?urlpath=lab/tree/notebooks/auto_examples/applications/plot_tomography_l1_reconstruction.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_tomography_l1_reconstruction.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_tomography_l1_reconstruction.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_tomography_l1_reconstruction.zip ` .. include:: plot_tomography_l1_reconstruction.recommendations .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_