Compressed Sensing Example

In this example we will go through a simple example from the field of Compressed Sensing. In addtion to RegularizedLeastSquares.jl, we will need the packages LinearOperatorCollection.jl, ImagePhantoms.jl, ImageGeoms.jl and Random.jl, as well as CairoMakie.jl for visualization. We can install them the same way we did RegularizedLeastSquares.jl.

Preparing the Inverse Problem

To get started, let us generate a simple phantom using the ImagePhantoms and ImageGeom packages:

using ImagePhantoms, ImageGeoms
N = 256
image = shepp_logan(N, SheppLoganToft())
size(image)

(256, 256)

This produces a 256x256 image of a Shepp-Logan phantom.

In this example, we consider a problem in which we randomly sample a third of the pixels in the image. Such a problem and the corresponding measurement can be constructed with the packages LinearOperatorCollection and Random:

We first randomly shuffle the indices of the image and then select the first third of the indices to sample.

using Random, LinearOperatorCollection
randomIndices = shuffle(eachindex(image))
sampledIndices = sort(randomIndices[1:div(end, 3)])

21845-element Vector{Int64}:
     2
     3
     4
     5
     6
     7
    13
    19
    22
    27
     ⋮
 65510
 65513
 65514
 65516
 65519
 65523
 65524
 65530
 65531

Afterwards we build a sampling operator which samples the image at the selected indices.

A = SamplingOp(eltype(image), pattern = sampledIndices , shape = size(image));

Then we apply the sampling operator to the vectorized image to obtain the sampled measurement vector

b = A*vec(image);

To visualize our image we can use CairoMakie:

using CairoMakie
function plot_image(figPos, img; title = "", width = 150, height = 150)
  ax = CairoMakie.Axis(figPos; yreversed=true, title, width, height)
  hidedecorations!(ax)
  heatmap!(ax, img)
end
fig = Figure()
plot_image(fig[1,1], image, title = "Image")
samplingMask = fill(false, size(image))
samplingMask[sampledIndices] .= true
plot_image(fig[1,2], image .* samplingMask, title = "Sampled Image")
resize_to_layout!(fig)
fig

As we can see in the right image, only a third of the pixels are sampled. The goal of the inverse problem is to recover the original image from this measurement vector.

Solving the Inverse Problem

To recover the image from the measurement vector, we solve the TV-regularized least squares problem:

\[\begin{equation} \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \lambda\vert\vert\mathbf{x}\vert\vert_{\text{TV}} . \end{equation}\]

For this purpose we build a TV regularizer with regularization parameter $λ=0.01$

using RegularizedLeastSquares
reg = TVRegularization(0.01; shape=size(image));

We will use the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) to solve our inverse problem. Thus, we build the corresponding solver

solver = createLinearSolver(FISTA, A; reg=reg, iterations=20);

and apply it to our measurement vector

img_approx = solve!(solver,b)

65536-element Vector{Float32}:
 7.5375277f-19
 1.4254111f-18
 3.4523988f-18
 8.367589f-18
 1.9255124f-17
 4.2910904f-17
 9.2002786f-17
 1.8753667f-16
 3.830775f-16
 8.49308f-16
 ⋮
 2.435243f-16
 1.2161483f-16
 5.90744f-17
 2.7788742f-17
 1.2509881f-17
 5.4832657f-18
 2.4100902f-18
 1.1232092f-18
 6.593036f-19

To visualize the reconstructed image, we need to reshape the result vector to the correct shape. Afterwards we can use CairoMakie again:

img_approx = reshape(img_approx,size(image));
plot_image(fig[1,3], img_approx, title = "Reconstructed Image")
resize_to_layout!(fig)
fig

This page was generated using Literate.jl.