Compressed Sensing in Julia

Yusri D. Heryanto

In this article, we will explore the application of compressed sensing using Julia. This post is inspired by a blog from Robert Taylor about compressed sensing in Python. Compressed sensing(CS) is a novel sensing paradigm that asserts one can recover signals or images from far fewer samples or measurements than conventional methods use.

CS relies on the principle of sparsity. Most natural signals, such as images and audio, are sparse when the signal is written on an appropriate basis. The commonly used bases for compression in images and audio are the Fourier or wavelet transform.

Mathematically, if a signal inline_formula not implemented is sparse in inline_formula not implemented basis, then instead of measuring y directly (inline_formula not implemented measurement). We can measure a randomly chosen measurement inline_formula not implementedwhere inline_formula not implemented can be dramatically fewer than inline_formula not implemented. The measurements inline_formula not implemented are given by

The matrix inline_formula not implemented is a sampling matrix. As mentioned above, we assumed inline_formula not implemented is sparse on some appropriate basis inline_formula not implemented:

Here, the variable inline_formula not implemented is a sparse vector. If the inline_formula not implemented basis is the Fourier basis, inline_formula not implemented can be seen as the inverse Fourier transform of s.

Then, the objective of compressed sensing is to find the sparsest vector inline_formula not implemented from the equation above. This is the optimization problem:

However, it is non-convex. Fortunately, we can relax the inline_formula not implemented problem to the inline_formula not implemented optimization problem:

The interested reader can read many papers and articles online such as "A User’s Guide to Compressed Sensing for Communications Systems".

One Dimensional Compressed sensing

First, we explore 1-dimensional compressed sensing. We use the example from Brunton and Kutz. The signal is a mixture of two signals with 2 different frequencies. For the sparse basis choice, we use Discrete Cosine Transform basis.

pkg"add Images FFTW StructuredOptimization Distributions Plots"

using Images

using FFTW

using Plots

using LinearAlgebra

using StructuredOptimization

using Distributions

# Create a signal which composed from 2 frequencies

n = 4096 # the number of time points

sampling_rate= n

t = range(0.0, step=1/sampling_rate, length=n)

y = sin.(2* 97 * pi * t) + sin.(2* 777 * pi * t) #freq 194 and 1554

freq = fftfreq(n,sampling_rate);

1+1

begin

	p1_ts =plot(t,y,legend=false,title = "Time domain signal")

	p2_ts = plot(t,y,xlim=(0.2,0.3),title = "Inset at [0.2,0.3]")

	plot(p1_ts,p2_ts,legend=false)

end

plot(freq,dct(y),xlim=(0,2000),legend=false,title="Frequency domain signal")

We showed you the signal in the time domain and in the frequency domain. Even though the signal in the time domain is not sparse but the signal in the frequency domain is sparse. This is called the incoherence. The incoherence of the signal is necessary conditions for signal recovery using CS.

# sample a small number of signal

n_sample= div(n,10)

idx = sort(rand(1:n,n_sample)) # sorting is necessary only for nice plotting

#the y_sample here is refer to x in equation above#

y_sample = y[idx];

plot(t[idx],y_sample,seriestype=:scatter,title="Sampling",label="Sampling point",legend=false)

plot!(t[idx],y_sample,seriestype=:line)

We sample only 10.0 percent of the signal.

Now, using CS theory, we can solve this problem using optimization:

However, we approached this problem using least absolute shrinkage and selection operator(LASSO) Algorithm:

s = Variable(n)

~s .= 0.0

lambda = 1e-1

# the function idct(s)[idx] here is same as  performing θ×s

@minimize ls(idct(s)[idx]-y_sample)+lambda*norm(s,1)

m_color = palette(:default)[1] #color for recovered signal

t_color = palette(:default)[2] # color for complete signal

plot(plot(freq,dct(y),color=t_color,title="Complete signal (frequency domain)"),

	plot(freq,~s,color=m_color,title="LASSO approach"),

	xlim=(0,2000),link=:both,legend=false)

plot(plot(t,y,color=t_color,title="Complete signal(time domain) "),

	plot(t,idct(~s,1),color=m_color,title="LASSO approach"),

	legend=false)

plot(plot(t,y,title="Complete signal inset",color=t_color),

	plot(t,idct(~s,1),title="LASSO approach inset",color=m_color),

	xlim=(0.2,0.23),legend=false)

The result is pretty good for the 10.0 percent of data. We can see that the recovery on the boundary is not quite good because violation of the periodic boundary condition requirements of the cosine transform.

Image reconstruction using compressed sensing

Image reconstruction can be seen as 2D or 3D reconstruction. Our example is using an image with three RGB channel colors. We have 2 approaches :(1) First, we can perform the 2D-DCT on each channel or we can perform 3D-DCT. We will show you both approach.

img = load(
bedroom_in_arles.jpg)[1:512,1:512] # load image with size 512*512*3

# We create this helper function for random sampling without replacement

function rand_wo_replacement(iterable,n_sample)

	res=Array{eltype(iterable),1}(undef,n_sample)

	i= 1 

	while i <= n_sample

		x = rand(iterable)

		if x ∈ res

			continue

		else 

			res[i] = x

		end

		i += 1

	end

	return res

end;

ny,nx=size(img)

len_img = length(img)

# we will sample about 5% of the image

n_sample_img= round(Int,len_img*0.05);

sample_idx=rand_wo_replacement(1:len_img,n_sample_img)

# We create a mask matrix

mask = zeros(size(img))

mask[sample_idx] .= 1.0

sample_img= mask .* img;

[img Gray.(mask) sample_img]

Here we have the original image, the mask, and the sampled image. We sample about 5.0 percent of the original image.

For the first approach, we extract each channel of the sampled image and flatten it. Then we perform compressed sensing on each channel.

2D Compressed sensing on each channel

#sample the image and flatten it

sample_img_array=img[sample_idx];

# Extract each channel of the RGB matrix. This will produced 3*length_of_image matrix. Each row represents each channel.

chanview = convert(Matrix{Float64},channelview(sample_img_array));

r_sample_img_arr= chanview[1,:] # R channel of the sampled image

g_sample_img_arr= chanview[2,:] # G channel of the sampled image

b_sample_img_arr= chanview[3,:]; # B channel of the sampled image

begin

	lambda_2= 3e-2

	r_vx2 = Variable(ny,nx) 

	g_vx2 = Variable(ny,nx) 

	b_vx2 = Variable(ny,nx) 

	@minimize ls(idct(r_vx2)[sample_idx]-r_sample_img_arr)+lambda_2*norm(r_vx2,1) with PANOC()

	@minimize ls(idct(g_vx2)[sample_idx]-g_sample_img_arr)+lambda_2*norm(g_vx2,1) with PANOC()

	@minimize ls(idct(b_vx2)[sample_idx]-b_sample_img_arr)+lambda_2*norm(b_vx2,1) with PANOC()

end

To see the recovered signal, we inverse transform our variable(i.e. idct(~r_vx2),idct(~g_vx2),idct(~b_vx2))

[colorview(RGB,idct(~r_vx2),idct(~g_vx2),idct(~b_vx2)) img]

3D Compressed sensing on RGB image

For the second approach, we perform 3 Dimensional DCT transform on the 3D array of the image

# Create 3D array from the sampled image

mat = permutedims(channelview(sample_img),(2,3,1));

# Find the linear index of the 3D array where is not zero

sample_idx3=LinearIndices(mat)[findall(x-> x>0, mat)]

# flatten the 3D array

sample_img_arr = mat[sample_idx3];

begin

	lambda_3= 3e-3

	vx3 = Variable(ny,nx,3)

	@minimize ls(idct(vx3)[sample_idx3]-sample_img_arr)+lambda_3*norm(vx3,1) with PANOC()

end

The get the image back from vx3, we inverse transform vx3 into n x n x 3 array and then permute the dimensions of array A to show as RGB image.

[colorview(RGB,permutedims(idct(~vx3),(3,1,2))) img]

Here we performed 2 approaches for decompressed sensing of RGB image. The result is recognizable. It shows you the capability of compressed sensing to recover the signal from far less sample.

Reference

[1] Brunton, Steven L., and Jose N. Kutz. Data-driven science and engineering: machine learning, dynamical systems, and control. Cambridge, United Kingdom New York, NY: Cambridge University Press, 2019. Print.

[2] Taylor, R. (2016, July 20). Compressed Sensing in Python [Web log post]. Retrieved October 12, 2020, from http://www.pyrunner.com/weblog/2016/05/26/compressed-sensing-python/

[3] Hayashi, K., Nagahara, M., & Tanaka, T. (2013). A User’s Guide to Compressed Sensing for Communications Systems. IEICE Transactions on Communications, E96.B(3), 685–712. https://doi.org/10.1587/transcom.e96.b.685