Lensless imaging with the Raspberry Pi and Python

7 min readDec 6, 2021

Now that we have built a lensless camera and taken some measurements, let’s see how to reconstruct an image. This task corresponds to what is typically referred to as an inverse problem. There are several approaches to solving inverse problems; in this tutorial we briefly present ADMM (Alternating Direction Method of Multipliers). Moreover, you will be asked to implement Accelerated Proximal Gradient Descent (APGD) and Primal Dual Splitting (PDS) as well as several different linear operators using Pycsou — a specialized package for solving linear inverse problems.

Source: https://waller-lab.github.io/DiffuserCam/tutorial/building_guide.pdf

This tutorial assumes that you have followed this one, and have downloaded and installed the LenslessPiCam repository.

Note that we have denoted exercises for students as such:

This is an exercise.

If you are an instructor / trying to replicate this tutorial, feel free to send an email to eric[dot]bezzam[at]epfl[dot]ch.

1) Introduction to lensless cameras

Conventional (top) vs. lensless (bottom) imaging. Image by author, inspired by Fig 4. of this paper.

Lenses are a key component of conventional cameras as they efficiently focus light onto an image plane. However, given a digital sensor the lens can be replaced with a different encoding element and a computational algorithm (see image above).

There are a couple of benefits of lensless cameras:

Thin / flexible form factor, as a potentially thick lens and tube configuration can be replaced by a thin mask (of arbitrary shape) closely placed to the sensor.
3-D imaging and computational refocusing, as we obtain an unfocused image from a lensless camera.

For an overview of lensless cameras, their benefits, and shortcoming, we refer to this overview paper.

2) Formulating the inverse problem

A key component of lensless cameras is formulating the inverse problem for recovering an image from our measurement. In the original DiffuserCam tutorial, they present and provide implementations for two inverse problem solvers: gradient descent and ADMM. We recommend reading their algorithms guide if you would like some insight behind setting up our inverse problem and how they solve it via gradient descent and ADMM. We briefly describe the main idea below.

Identifying the forward model

We will be making the same assumptions for the forward model (Section 2.1), namely that we have a linearly shift-invariant (LSI) system.

Linear shift-invariant (LSI) system. Source: https://waller-lab.github.io/DiffuserCam/tutorial/algorithm_guide.pdf

This LSI assumption is what allows us to measure a single PSF, as our measurement will be a sum of weighted and shifted versions of our PSF. Furthermore, our measurement can be modeled as a convolution between the PSF and the object in the scene.

An important component for ensuring LSI is the cropping of the PSF we did when constructing the camera. If we didn’t have this cropping, new caustic patterns would appear if we were to shift the PSF, and we could no longer model our measurement as superpositions of the same PSF weighted and shifted.

Therefore, our forward model consists of a convolution followed by a cropping.

We can express these operations as matrix-vector products by stacking the columns of the measurement y and the scene x. Moreover, applying C and P correspond to applying cropping and convolution with our PSF, which we combine into a single operation H (this is the operation we estimate with our PSF measurement).

Inverting the problem

We measure y but we would like to recover x. A common way of obtaining x would be by inverting H:

However, H in our case is not invertible. One reason being that C (cropping) results in a loss of information.

For such scenarios, we can resort to iterative algorithms that minimize a data fidelity term.

Various priors and regularization terms can be incorporated in order to speed-up / guide the convergence of the iterative algorithm. For example, we can safely assume that our recovered image must be non-negative:

In the original tutorial, they apply gradient descent and ADMM (Alternating Direction Method of Multipliers) to minimize this term with the non-negativity constraint. We refer to the linked notebooks, as they provide an in-depth explanation of the two approaches along with code.

In the LenslessPiCam library, we have implemented an object-oriented version of ADMM, which we have adapted from the original tutorial so that it can reconstruct RGB images. We have also done some optimization, e.g. pre-computing some matrices and using the real-valued FFT. You can run the algorithm like so:

python scripts/recon/admm.py

File paths and default parameters are set in configs/admm_thumbs_up.yaml.

We can use this approach as a baseline for the other techniques we will implement with Pycsou.

3) Solving inverse problems with Pycsou (exercises)

For the remainder of the tutorial, we focus on using Pycsou, which provides a modular interface for defining inverse problems and contains several already implemented solvers and penalties.

The first step is to install Pycsou inside our virtual environment (if not done already). This can be done with pip.

(lensless_env) pip install pycsou

The first approach we will implement is accelerated proximal gradient descent (APGD).

Exercise: implement a penalized least-square problem with Tikhonov regularization (i.e. ridge regression):

where y is the (flattened) grayscale raw measurements obtained with LenslessPiCam, H is a discrete convolution operator modeling the blur introduced by LenslessPiCam’s PSF, x is the (flattened) grayscale image to recover, and lambda is a penalty parameter. Play with various penalty parameters. You may find the following Pycsou classes useful in your implementation: pycsou.opt.proxalgs.APGD, pycsou.func.penalty.SquaredL2Norm, and pycsou.linop.conv.Convolve2D. You can also use this template for loading the PSF and raw data. Be sure to use the --gray flag for grayscale data. See also Pycsou's documentation for an example describing how to use APGD.

We now modify the regularization term.

Exercise: perform the grayscale image reconstruction with a LASSO problem, i.e. replacing the squared L2 norm in the penalty term of (1) by an L1 norm. Again, investigate different penalty parameters. You may find the following Pycsou class useful in your implementation: pycsou.func.penalty.L1Norm.

We now try the same prior as in the original DiffuserCam tutorial.

Exercise: perform the grayscale image reconstruction with a non-negative least-square problem, i.e. replacing the squared L2 norm in the penalty term of (1) by a non-negativity prior. You may find the following Pycsou class useful in your implementation: pycsou.func.penalty.NonNegativeOrthant.

We now impose sparsity within a specific domain by applying a linear transformation.

Exercise: perform the grayscale image reconstruction with a generalized LASSO problem with the Type 2 Discrete Cosine Transform (DCT) as a sparsifying transform, i.e. replacing the square L2 norm in the penalty term of (1) by an L1 norm composed with the 2D Type 2 DCT transform. Pycsou does not support yet the DCT transform, so you will have to implement it yourself by sub-classing LinearOperator (the base class for linear operators in Pycsou) and implement the __call__ and adjoint methods (see Pycsou's documentation for an example). You may find the following SciPy classes useful in your implementation: scipy.fft.dctn and scipy.fft.idctn. Note that with the right normalization mode, the DCT of Type 2 is an orthonormal operator (see the documentation notes of scipy.fft.dct for more information on the topic). Use this fact to propose a change of variable that transforms the generalized LASSO problem with penalty |DCT_2(x)|_1 into a regular LASSO problem. Solve the LASSO problem resulting from this change of variable via pycsou.opt.proxalgs.APGD.

We now combine multiple penalty terms, first with the primal-dual splitting method (PDS), of which ADMM is a special case.

Exercise: perform the grayscale image reconstruction with a penalized least-square problem with composite penalty term involving a total variation (TV) and non-negativity prior. This means replacing the squared L2 norm in the penalty term of (1) by the sum of a non-negativity prior and a weighted L1 norm composed with the finite-difference 2D gradient operator. You may find the following Pycsou classes useful in your implementation: pycsou.linop.diff.Gradient and pycsou.opt.proxalgs.PDS (see Pycsou’s documentation for an example). Note that PDS is much slower than APGD, and may require at least 5000 iterations to converge.

To accelerate the previous reconstruction, we replace the L1 norm in the TV penalty by a smoothed L1 norm known as the Huber loss.

Exercise: Implement the Huber norm by sub-classing DifferentialFunctional and providing the __call__ and jacobianT methods. Do not forget to provide the Lipschitz constant of the derivative via the diff_lipschitz_cst argument of the DifferentiableFunctional.__init__ method. See Pycsou’s documentation for an example.

Now that we have implemented a whole slew of techniques in grayscale, let’s extend them to RGB!

Exercise: extend your implementations to RGB. You can take inspiration from ADMM which updates all channels simultaneously using array operations.

Finally, we can significantly speed up the Pycsou implementations by taking advantage of the fact that our forward operator (PSF) is fixed and that we are working with image intensities. This requires replacing pycsou.linop.conv.Convolve2D with a custom operator.