{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Compressive sensing: tomography reconstruction with L1 prior (Lasso)\n\nThis example shows the reconstruction of an image from a set of parallel\nprojections, acquired along different angles. Such a dataset is acquired in\n**computed tomography** (CT).\n\nWithout any prior information on the sample, the number of projections\nrequired to reconstruct the image is of the order of the linear size\n``l`` of the image (in pixels). For simplicity we consider here a sparse\nimage, where only pixels on the boundary of objects have a non-zero\nvalue. Such data could correspond for example to a cellular material.\nNote however that most images are sparse in a different basis, such as\nthe Haar wavelets. Only ``l/7`` projections are acquired, therefore it is\nnecessary to use prior information available on the sample (its\nsparsity): this is an example of **compressive sensing**.\n\nThe tomography projection operation is a linear transformation. In\naddition to the data-fidelity term corresponding to a linear regression,\nwe penalize the L1 norm of the image to account for its sparsity. The\nresulting optimization problem is called the `lasso`. We use the\nclass :class:`~sklearn.linear_model.Lasso`, that uses the coordinate descent\nalgorithm. Importantly, this implementation is more computationally efficient\non a sparse matrix, than the projection operator used here.\n\nThe reconstruction with L1 penalization gives a result with zero error\n(all pixels are successfully labeled with 0 or 1), even if noise was\nadded to the projections. In comparison, an L2 penalization\n(:class:`~sklearn.linear_model.Ridge`) produces a large number of labeling\nerrors for the pixels. Important artifacts are observed on the\nreconstructed image, contrary to the L1 penalization. Note in particular\nthe circular artifact separating the pixels in the corners, that have\ncontributed to fewer projections than the central disk.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The scikit-learn developers\n# SPDX-License-Identifier: BSD-3-Clause\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy import ndimage, sparse\n\nfrom sklearn.linear_model import Lasso, Ridge\n\n\ndef _weights(x, dx=1, orig=0):\n    x = np.ravel(x)\n    floor_x = np.floor((x - orig) / dx).astype(np.int64)\n    alpha = (x - orig - floor_x * dx) / dx\n    return np.hstack((floor_x, floor_x + 1)), np.hstack((1 - alpha, alpha))\n\n\ndef _generate_center_coordinates(l_x):\n    X, Y = np.mgrid[:l_x, :l_x].astype(np.float64)\n    center = l_x / 2.0\n    X += 0.5 - center\n    Y += 0.5 - center\n    return X, Y\n\n\ndef build_projection_operator(l_x, n_dir):\n    \"\"\"Compute the tomography design matrix.\n\n    Parameters\n    ----------\n\n    l_x : int\n        linear size of image array\n\n    n_dir : int\n        number of angles at which projections are acquired.\n\n    Returns\n    -------\n    p : sparse matrix of shape (n_dir l_x, l_x**2)\n    \"\"\"\n    X, Y = _generate_center_coordinates(l_x)\n    angles = np.linspace(0, np.pi, n_dir, endpoint=False)\n    data_inds, weights, camera_inds = [], [], []\n    data_unravel_indices = np.arange(l_x**2)\n    data_unravel_indices = np.hstack((data_unravel_indices, data_unravel_indices))\n    for i, angle in enumerate(angles):\n        Xrot = np.cos(angle) * X - np.sin(angle) * Y\n        inds, w = _weights(Xrot, dx=1, orig=X.min())\n        mask = np.logical_and(inds >= 0, inds < l_x)\n        weights += list(w[mask])\n        camera_inds += list(inds[mask] + i * l_x)\n        data_inds += list(data_unravel_indices[mask])\n    proj_operator = sparse.coo_matrix((weights, (camera_inds, data_inds)))\n    return proj_operator\n\n\ndef generate_synthetic_data():\n    \"\"\"Synthetic binary data\"\"\"\n    rs = np.random.RandomState(0)\n    n_pts = 36\n    x, y = np.ogrid[0:l, 0:l]\n    mask_outer = (x - l / 2.0) ** 2 + (y - l / 2.0) ** 2 < (l / 2.0) ** 2\n    mask = np.zeros((l, l))\n    points = l * rs.rand(2, n_pts)\n    mask[(points[0]).astype(int), (points[1]).astype(int)] = 1\n    mask = ndimage.gaussian_filter(mask, sigma=l / n_pts)\n    res = np.logical_and(mask > mask.mean(), mask_outer)\n    return np.logical_xor(res, ndimage.binary_erosion(res))\n\n\n# Generate synthetic images, and projections\nl = 128\nproj_operator = build_projection_operator(l, l // 7)\ndata = generate_synthetic_data()\nproj = proj_operator @ data.ravel()[:, np.newaxis]\nproj += 0.15 * np.random.randn(*proj.shape)\n\n# Reconstruction with L2 (Ridge) penalization\nrgr_ridge = Ridge(alpha=0.2)\nrgr_ridge.fit(proj_operator, proj.ravel())\nrec_l2 = rgr_ridge.coef_.reshape(l, l)\n\n# Reconstruction with L1 (Lasso) penalization\n# the best value of alpha was determined using cross validation\n# with LassoCV\nrgr_lasso = Lasso(alpha=0.001)\nrgr_lasso.fit(proj_operator, proj.ravel())\nrec_l1 = rgr_lasso.coef_.reshape(l, l)\n\nplt.figure(figsize=(8, 3.3))\nplt.subplot(131)\nplt.imshow(data, cmap=plt.cm.gray, interpolation=\"nearest\")\nplt.axis(\"off\")\nplt.title(\"original image\")\nplt.subplot(132)\nplt.imshow(rec_l2, cmap=plt.cm.gray, interpolation=\"nearest\")\nplt.title(\"L2 penalization\")\nplt.axis(\"off\")\nplt.subplot(133)\nplt.imshow(rec_l1, cmap=plt.cm.gray, interpolation=\"nearest\")\nplt.title(\"L1 penalization\")\nplt.axis(\"off\")\n\nplt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)\n\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.21"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}