Reconstruction Methods in THz Single-pixel Imaging

TL;DR

This paper explores advanced image reconstruction techniques in THz single-pixel imaging, comparing regularization methods, addressing nonlinear calibration and phase recovery challenges, and discussing future research directions.

Contribution

It provides a comprehensive comparison of regularization techniques and discusses nonlinear reconstruction issues specific to THz single-pixel cameras.

Findings

Regularization techniques significantly affect reconstruction quality.

Calibration and phase recovery pose nonlinear challenges.

Future research directions are outlined for improved imaging.

Abstract

The aim of this paper is to discuss some advanced aspects of image reconstruction in single-pixel cameras, focusing in particular on detectors in the THz regime. We discuss the reconstruction problem from a computational imaging perspective and provide a comparison of the effects of several state-of-the art regularization techniques. Moreover, we focus on some advanced aspects arising in practice with THz cameras, which lead to nonlinear reconstruction problems: the calibration of the beam reminiscent of the Retinex problem in imaging and phase recovery problems. Finally we provide an outlook to future challenges in the area.

Figures

Equations

$$\mathbf{y}=\mathbf{Ap},$$

$$\Vert \mathbf{\alpha }{\Vert }_{{\ell }_0}=|\{k:{\alpha }_k\ne 0\}|,$$

$$|\{k:|{\alpha }_k|\ge ϵ\}|,$$

$$\underset{\mathbf{\alpha }}{\min }\Vert \mathbf{\alpha }{\Vert }_{{\ell }_1}\text{s.t.}\mathbf{y}=\mathbf{A\Psi }\mathbf{\alpha }$$

$$\underset{\mathbf{\alpha }}{\min }\Vert \mathbf{\alpha }{\Vert }_{{\ell }_1}\text{s.t.}\mathbf{y}=\mathbf{A\Psi }\mathbf{\alpha }\text{and}\mathbf{\Psi }\mathbf{\alpha }\ge 0.$$

$$\begin{array}{rl}\Vert \mathbf{p}{\Vert }_{{\ell }_1}-\Vert {\mathbf{p}}^0{\Vert }_{{\ell }_1}& =⟨1,\mathbf{p}-{\mathbf{p}}^0⟩=⟨\mathbf{t},\mathbf{A}(\mathbf{p}-{\mathbf{p}}^0)⟩\\ & =⟨\mathbf{t},\mathbf{Ap}-\mathbf{y}⟩\le \Vert \mathbf{t}{\Vert }_{{\ell }_2}\Vert \mathbf{Ap}-\mathbf{y}{\Vert }_{{\ell }_2}\end{array}$$

$$\underset{\mathbf{p}\ge 0}{\min }\Vert \mathbf{y}-\mathbf{Ap}{\Vert }_{{\ell }_2}^2.$$

$$\underset{\mathbf{\alpha }}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{A\Psi }\mathbf{\alpha }{\Vert }_{{\ell }_2}^2+\lambda \Vert \mathbf{\alpha }{\Vert }_{{\ell }_1}\text{s.t.}\mathbf{\Psi }\mathbf{\alpha }\ge 0,$$

$$\underset{\mathbf{p}}{\min }\Vert \nabla \mathbf{p}{\Vert }_{{\ell }_1}\text{s.t.}\mathbf{y}=\mathbf{Ap}$$

$$\underset{\mathbf{p}}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{Ap}{\Vert }_{{\ell }_2}^2+\lambda \Vert \nabla \mathbf{p}{\Vert }_{{\ell }_1},$$

$$\underset{\mathbf{p},\mathbf{w}}{\min }\Vert \nabla \mathbf{p}-\mathbf{w}{\Vert }_{{\ell }_1}+\beta \Vert \nabla \mathbf{w}{\Vert }_{{\ell }_1}\text{s.t.}\mathbf{y}=\mathbf{Ap},$$

$$\underset{\mathbf{p}}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{Ap}{\Vert }_{{\ell }_2}^2+\lambda R(\mathbf{\Phi }\mathbf{p}),$$

$$\underset{\mathbf{p}}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{Ap}{\Vert }_{{\ell }_2}^2+\tilde{\lambda }(R(\mathbf{\Phi }\mathbf{p})-{\mathbf{q}}^k\cdot \mathbf{p})$$

$${\mathbf{p}}_{k+1}\in \text{arg}\underset{\mathbf{p}}{\min }\frac{\mu }{2}\Vert \mathbf{y}-\mathbf{Ap}{\Vert }_{{\ell }_2}^2+R(\mathbf{\Phi }\mathbf{p})+{\mathbf{r}}^k\cdot (\mathbf{y}-\mathbf{Ap})$$

$${\mathbf{r}}^{k+1}={\mathbf{r}}^k+\mathbf{y}-\mathbf{A}{\mathbf{p}}^{k+1}.$$

$${\text{prox}}_{{\ell }_1}(f)=\text{arg}\underset{x}{\min }\frac{1}{2}\Vert x-f{\Vert }_{{\ell }_2}^2+\Vert x{\Vert }_{{\ell }_1}$$

$$\underset{\mathbf{p}}{\min }F(\mathbf{p})+G(\mathbf{L}\mathbf{p}),$$

$$\underset{\mathbf{p}}{\min }\underset{\mathbf{q}}{\max }F(\mathbf{p})+⟨\mathbf{L}\mathbf{p},\mathbf{q}⟩-G^{\ast }(\mathbf{q}),$$

$$\mathbf{p}=\mathbf{d}\odot \mathbf{x},$$

$$\underset{\mathbf{x}\in [0,1{]}^n,\mathbf{d}\ge 0}{\min }\Vert \mathbf{y}-\mathbf{A}(\mathbf{d}\odot \mathbf{x}){\Vert }_{{\ell }_2}^2+\lambda \cdot r(\mathbf{x},\mathbf{d})$$

$$\mathbf{d}\odot \mathbf{x}=\frac{1}{\sqrt{n}}{\mathbf{F}}^{-1}[(\mathbf{Fd})\ast (\mathbf{Fx})]=\frac{1}{\sqrt{n}}{\mathbf{F}}^{-1}[\mathbf{l}\ast \mathbf{r}]$$

$$\underset{\mathbf{x}\in [0,1{]}^n,\mathbf{d}\ge 0}{\min }\Vert \mathbf{p}-\mathbf{d}\odot \mathbf{x}{\Vert }_{{\ell }_2}^2+\lambda \cdot r(\mathbf{x},\mathbf{d}).$$

$$\underset{\mathbf{x}\in [0,1{]}^n,\mathbf{g}>ϵ}{\min }\Vert \mathbf{p}-\mathbf{g}\odot \mathbf{Ax}{\Vert }_{{\ell }_2}^2+\lambda \cdot r(\mathbf{x},\mathbf{g}).$$

$$\underset{\mathbf{l}\ge \mathbf{s}}{\min }\frac{1}{2}\Vert \mathbf{l}-\mathbf{s}{\Vert }_{{\ell }_2}^2+\frac{\alpha }{2}\Vert \nabla (\mathbf{l}-\mathbf{s}){\Vert }_{{\ell }_2}^2+\frac{\beta }{2}\Vert \nabla \mathbf{l}{\Vert }_{{\ell }_2}^2.$$

$$\underset{\mathbf{l}\ge \mathbf{s},\mathbf{r}\le 0}{\min }\frac{1}{2}\Vert \mathbf{l}-\mathbf{s}+\mathbf{r}{\Vert }_{{\ell }_2}^2+\alpha \Vert \nabla \mathbf{r}{\Vert }_{{\ell }_1}+\frac{\beta }{2}\Vert \nabla \mathbf{l}{\Vert }_{{\ell }_2}^2.$$

$$\underset{\mathbf{l}\ge \mathbf{s},\mathbf{r}\le 0}{\min }\Vert \mathbf{y}-\mathbf{A}(e^{\mathbf{l}+\mathbf{r}}){\Vert }_{{\ell }_2}^2+\alpha \Vert \nabla \mathbf{r}{\Vert }_{{\ell }_1}+\frac{\beta }{2}\Vert \nabla \mathbf{l}{\Vert }_{{\ell }_2}^2.$$

$$\underset{\mathbf{r}\le 0,\mathbf{l}\in {\mathbb{R}}^n}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{A}{e}^{\mathbf{r}+\mathbf{l}}{\Vert }_{{\ell }_2}^2+\lambda \cdot \tilde{r}(\mathbf{r},\mathbf{l})$$

$${\mathbf{r}}^{k+1/2}$$

$${\mathbf{l}}^{k+1/2}$$

$$\frac{1}{2\tau }\Vert \mathbf{r}-{\mathbf{r}}^{k+1/2}{\Vert }_{{\ell }_2}^2+{\lambda }_r\Vert \nabla r{\Vert }_{{\ell }_1}$$

$$\frac{1}{2\tau }\Vert \mathbf{l}-{\mathbf{l}}^{k+1/2}{\Vert }_{{\ell }_2}^2+{\lambda }_l\Vert \nabla l{\Vert }_{{\ell }_2}^2.$$

Full text

Moreover, we focus on some advanced aspects arising in practice with THz cameras, which lead to nonlinear reconstruction problems: the calibration of the beam reminiscent of the Retinex problem in imaging and phase recovery problems. Finally we provide an outlook to future challenges in the area.

1 Introduction

Imaging science has been a strongly evolving field in the last century, with a lot of interesting developments concerning devices, measurement strategies and computational approaches to obtain high-quality images. A current focus concerns imaging from undersampled data in order to allow novel developments towards dynamic and hyperspectral imaging, where time restrictions forbid to acquire full samplings. In order to compensate for the sampling either physical models (e.g. concerning motion in dynamic imaging, cf. [BDS18]) or a-priori knowledge about the images to be reconstructed are used. In applications, where one has sufficient freedom to choose the undersampling patterns, the paradigm of compressed sensing is particularly popular. It is based on minimizing coherence between measurements (cf. [D*+*06]), often achieved by random sampling (cf. [Bar07]).

Single-pixel imaging, or more precisely single-detector imaging, is one of the most interesting developments in compressed sensing (cf. [DDT*+*08] [CCT*+*08] [EGP18] [SSPB18]). It is based on using a single detector with multiple random masks in order to achieve the desired resolution. The ideal model is to have the mask realize a grid on the detector region with subpixels of the desired image resolution, which are either open or closed with a certain probability. The detector integrates the light passing through the open subpixels in the mask. Note that the image at desired resolution might be acquired by scanning in a deterministic way, in the easiest setting with a single subpixel open at each shot. However, the light intensity obtained from a single subpixel might not be sufficient to obtain a reasonable signal-to-noise ratio and obviously the mechanical scanning times may strongly exceed the potential times of undersampling with masks having multiple open subpixels.

This idea is particularly relevant in applications where detectors are expensive or difficult to miniaturize such as imaging in the THz range (cf. [CCT*+*08] [WSM*+*14] [AHJH15] [AFJH17]), which is our main source of motivation. The random masks are achieved by a spatial light modulator, which may however deviate from the ideal setting in practical applications due to the following effects:

• •

Beam calibration: as in many practical imaging approaches, the lighting beam is not homogeneous and needs to be corrected, a problem reminiscent of the classical Retinex problem (cf. [LM71]).

• •

Diffraction: deviations between the object and image plane may cause the need to consider diffraction and take care of out-of phase effects, which complicates the inversion problem and

• •

*Motion: * The object to be imaged may move between consecutive shots, i.e. while changing the masks, which leads to well-known motion blur effects in reconstructions.

In the following we will discuss such issues in image reconstruction arising in THz single pixel cameras. We start with the basic modeling of the image formation excluding calibration and further effects in Section 2, where we also discuss several regularization models. In Section 3 we discuss some computational methods to efficiently solve the reconstruction problem and in particular compare the results of some state-of-the art regularization techniques. Then we discuss the challenges that arise when applying the single-pixel imaging approach in a practical setup with THz detectors and proceed to calibration models for beams in Section 5, where we provide a relation to the classical Retinex problem and discuss the particular complications arising due to the combination with the single pixel camera. In Section 6 we discuss the phase reconstruction problem and compare several state-of-the art approaches for such. Finally we conclude and discuss further challenges in Section 7.

2 Compressed Sensing and Reconstruction in Single Pixel Cameras

In the following we discuss some basic aspects of the image reconstruction in single pixel cameras. We are interested in reconstructing a two-dimensional image on the subpixel grid of size $d_{1}\times d_{2}$, i.e. $\mathbf{p}\in\mathbb{R}^{n}$ where $n=d_{1}\cdot d_{2}$. The simplest model of the image formation process is to have each measurement $y_{i}$, $i=1,\ldots,m$, as the sum of subpixel values $p_{j}$ for those $j$ corresponding to the open subpixels in the $i$-th mask. This means we can write

[TABLE]

with a matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ whose entries are only zeros and ones. Again the non-zero entries in each row of $\mathbf{A}$ correspond to the open subpixels of the respective mask (compare Figure 1 for a schematic overview of such a setup).

Choosing $m\geq n$ deterministic masks appropriately one could guarantee that $\mathbf{A}$ is invertible and simply solve the linear reconstruction problem. However, in practice one would like to reconstruct a reasonable image $\mathbf{p}$ from a number of measurements $m$ significantly smaller than $n$. Single pixel cameras are hence developed in the paradigm of compressed sensing and the natural way to realize appropriate measurements is to choose the masks randomly. In most systems such as the ones we consider here, each entry of $\mathbf{A}$ is chosen independently as a binary random variable with fixed expectation. Combining such approaches with appropriate prior information in terms of sparsity leads to compressed sensing schemes that can be analyzed rigorously (cf. [CP11a] [LLKR15]).

2.1 Compressed Sensing Techniques

A key motivation for compressed sensing comes from the fact that in many cases $d_{1}\times d_{2}$ images are compressible and can be (approximately) sparsely represented as $\mathbf{p}=\mathbf{\Psi}\boldsymbol{\alpha}$ with $\boldsymbol{\alpha}\in\mathbb{R}^{n}$ and a particular basis $\mathbf{\Psi}\in\mathbb{R}^{n\times n}$, e.g. wavelets, or overcomplete systems $\mathbf{\Psi}\in\mathbb{R}^{n\times N}$ with $N>n$ such as shearlets (cf. [GL07]). By this we mean that

[TABLE]

respectively

[TABLE]

for small $\epsilon>0$, is considerably smaller then the ambient dimension $n$. Ideally, one would thus solve the problem of minimizing $\|\boldsymbol{\alpha}\|_{\ell_{0}}$ subject to the linear constraint $\mathbf{A}\mathbf{\Psi}\boldsymbol{\alpha}=\mathbf{y}$, which is however an NP-hard combinatorial problem.

One of the fundamental results that initiated the field of compressed sensing (cf. [CRT06, Don06, EK12, FR13]) is that under additional assumptions on the map $\mathbf{A}\mathbf{\Psi}$ the convex relaxation

[TABLE]

recovers exactly (in the noiseless setting) the unknown $\boldsymbol{\alpha}$ yielding the correct image $\mathbf{p}$. A common additional assumption is that the image itself has non-negative pixel intensities, $\mathbf{p}\geq 0$ such that problem (2) is extended to:

[TABLE]

If $\mathbf{\Psi}=\text{Id}$ and if the row span of $\mathbf{A}$ intersects the positive orthant, meaning that there exists a vector $\mathbf{t}$ such that $\mathbf{A}^{*}\mathbf{t}>0$, the measurement matrix $\mathbf{A}$ itself already assures that the $\ell_{1}$-norm $\|\mathbf{p}\|_{\ell_{1}}$ of a feasible $\mathbf{p}=\boldsymbol{\alpha}$ in (3) is equal (or close in the noisy case) to $\|\mathbf{p}^{0}\|_{\ell_{1}}$ where $\mathbf{p}^{0}$ is the unknown image to recover. To see this, let us assume exemplary that we find a vector $\mathbf{t}$ such that $\mathbf{A}^{*}\mathbf{t}=\mathbf{1}$ is the all-one vector and $\mathbf{y}=\mathbf{A}\mathbf{p}^{0}$. Then

[TABLE]

and hence, it is enough to minimize the residual over $\mathbf{p}\geq 0$ and replace (3) by a simple non-negative least squares (NNLS) problem:

[TABLE]

Indeed, that under these assumptions $\ell_{1}$-minimization reduces to a feasibility problem has observed already in prior work [BEZ08, WXT11]. In particular the setting of random binary masks has been investigated in [SH13, KJ17] and a considerably (partially-) derandomized result based on orthogonal arrays is discussed in [JKM18]. Although this explains to some extend why non-negativity and NNLS are very useful in certain imaging problems this does not easily extends to generic dictionaries $\mathbf{\Psi}$.

Hence, coming back to (3), in the case of noisy data knowledge about the expected solution should be included. Hence one usually rather solves

[TABLE]

where $\lambda>0$ is an appropriate regularization parameter.

2.2 Total Variation Regularization and Related Methods

While simple $\ell_{1}$-regularization is a common approach used for the theory of compressed sensing, for image reconstruction this choice as data fitting term has some drawbacks in practice, e.g. in wavelets the reconstructions may suffer from artifacts due to rectangular structures in their constructions. In more advanced models like shearlets (cf. [KL12]) the visual quality of reconstructions is improved, but the computational overhead in such large dictionaries may become prohibitive. A much more popular approach in practical image reconstructions are total variation based methods such as minimizing

[TABLE]

or the penalty version

[TABLE]

potentially with additional non-negativity constraints. Total variation methods in compressed sensing have been investigated recently in further detail (cf. [NW13b] [NW13a] [Poo15] [CX15] [KKS17]).

As in wavelet systems, simple approaches in total variation may suffer from rectangular grid artefacts. This is the case in particular if the straight $\ell_{1}$-norm of $\nabla\mathbf{p}$ is used in the regularization (namely $\|\nabla\mathbf{p}\|_{\ell_{1,1}}$), which corresponds to an anisotropic total variation. It is well-known that such approaches promote rectangular structures aligned with the coordinate axis (cf. [CCN14]), but destroy round edges. An improved version is the isotropic version, which considers $\nabla\mathbf{p}\in\mathbb{R}^{n\times 2}$ (the rows corresponding to partial derivatives) and computes $\|\nabla\mathbf{p}\|_{\ell_{2,1}}$ as total variation. The isotropic total variation promotes round structures that are visually more appealing in natural images, which can be made precise in a continuum limit.

A remaining issue of total variation regularization in applications is the so-called stair-casing phenomenon, which means that piecewise constant structures are promoted too strongly and thus gradual changes in an image are rather approximated by piecewise constants in a stair-like fashion. In order to cure such issues infimal convolutions of total variation and higher order functionals are often considered, the most prominent one being the total generalized variation (TGV) model (cf. [BKP10])

[TABLE]

where $\beta>0$ is a parameter to be chosen appropriately. Again in practice suitable isotropic variants of the $\ell_{1}$-norm are used, and in most approaches only the symmetrized gradient of the vector field ${\bf w}$ is used instead of the full gradient.

In the penalized form for noisy data used in practice, the models above can be written in the unified form

[TABLE]

with $R$ being an appropriately chosen seminorm and a suitable matrix ${\bf\Phi}$, e.g. the gradient or ${\bf\Phi}={\bf\Psi}^{-1}$ in the case of using a basis, or ${\bf\Phi}={\bf\Psi}^{*}$ for analysis sparsity. A common issue in problems of this form that aim to reduce variance is an increase of bias, which e.g. leads to a loss of contrast and small structures in total variation regularization. In order to improve upon this issue iterative regularization (cf. [BO13] [BB18]) can be carried out, in the simplest case by the Bregman iteration, which computes each iteration $\mathbf{p}^{k+1}$ as the minimizer of

[TABLE]

with subgradient ${\bf q}^{k}\in\partial R({\bf\Phi}\mathbf{p}^{k})$. In this way in each iteration a suitable distance to the last iterate is penalized instead of the original functional $R$, which is effectively a distance to zero. The parameter $\tilde{\lambda}$ is chosen much larger than the optimal $\lambda$ in the above problem, roughly speaking when carrying out $K$ iterations we have $\tilde{\lambda}=K\lambda.$ Observing from the optimality condition that ${\bf q}^{k}=\mu\mathbf{A}^{*}{\bf r}^{k}$ with $\mu=\frac{1}{\tilde{\lambda}}$, the Bregman iteration can be interpreted equivalently as the augmented Lagrangian method for the constrained problem of minimizing $R({\bf\Phi}\mathbf{p})$ subject to $\mathbf{y}=\mathbf{A}\mathbf{p}$, i.e.

[TABLE]

and

[TABLE]

While single pixel cameras are naturally investigated from a compressed sensing point of view, let us comment on some aspects of the problem when viewed as an inverse problem as other image reconstruction tasks in tomography or deblurring. First of all, we can see some common issues in computational techniques, which are needed due to the indirect relation between the image and the measurements and the fact that the matrix $\mathbf{A}$ is related to the discretization of an integral operator. On the other hand there is a significant difference between the single pixel setup and other image reconstruction problems in the sense that the adjoint operator is not smoothing, i.e. $\mathbf{A}^{*}$ has no particular structure. Thus, the typical source condition $\mathbf{A}^{*}{\bf{w}}\in{\bf{\Phi}}^{*}\partial R({\bf{\Phi}}\mathbf{p})$ for some ${\bf{w}}$ that holds for solutions of the variational problems (cf. [BB18]) does not imply smoothness of the subgradients as in other inverse problems.

3 Computational Image Reconstruction

In the following we discuss the effects of different regularization models on the image reconstruction quality. In order to efficiently compute solutions of the arising variational problems at reasonable image resolution, appropriate schemes to handle the convex but nondifferentiable terms are needed. It has become a standard approach in computational imaging of such problems to employ first-order splitting methods based on proximal maps that can be computed exactly, we refer to [BSS16] for an overview. The key idea is to isolate matrix-vector multiplications and local nonlinearities, for which the proximal map can be computed exactly or numerically efficient. A standard example is the $\ell_{1}$-norm, whose $\ell_{2}$-proximal map

[TABLE]

is given by soft shrinkage. A very popular approach in current computational imaging are first-order primal dual methods (cf. [CP11b] [ZBO11]) for computing minimizers of problems of the form

[TABLE]

with convex functionals $F$ and $G$ and a linear operator ${\bf L}$. The primal dual approach reformulates the minimization as a saddle point problem

[TABLE]

with $G^{*}$ being the convex conjugate of $G$. Then one iterates primal minimization with respect to $\mathbf{p}$ and dual maximization with respect to ${\bf q}$ with additional damping and possible extrapolation steps. Here we use a toolbox by Dirks [Dir16] for rapid prototyping with such methods to implement the models of the previous section.

We will use synthetic data to investigate the regularization methods described above. With this we can fully explore and understand the effects of named methods for different types of images. Here we consider the multiple regularization operators, namely Tikhonov, $\ell_{1}$- regularization, total variation (TV) and total generalized variation (TGV) [BKP10]. Moreover we compare the approach with a simple non-negative least squares approach that only enforces non-negativity of the image. Each dataset is composed of a ground truth structural image and a Gaussian kernel beam that is used as an illumination source. In this work we will explore two different ground truth images, that are illuminated by the shown beam. Structural images, beam and illuminated structures are shown in Figure 2.

An overview of basic reconstructions in the single detector framework using random binary masks is shown in Figure 3. Here we use two phantoms and two different sampling ratios. We display reconstructions with regularization parameters optimized for SSIM and observe the impact of the regularization for both choices of the sampling rate. In particular we see that for total variation type regularizations there is a strong improvement in visual image quality and that there is hardly any loss when going from full sampling to undersampling.

4 Challenges in Practical THz Single-Pixel Imaging

If compared to the visible region of the electromagnetic spectrum (VIS) the THz region is quite demanding from a physical point of view. Due to the large wavelength that is 2-3 orders of magnitude larger than for VIS radiation, the THz region is plagued by coherence effects and diffraction. This has to be mitigated already on the physical level by using specifically designed optical elements and techniques. As a rule of thumb, techniques and methods from the VIS region can be used but they need to be adapted. Also, lenses and mirrors are used but they are made of different materials or they tend to be bigger. For this reason the THz region is also sometimes called the quasi-optical or Gaussian region. It is called the Gaussian region because beams in the THz region almost always are of Gaussian shape. This, in turn, means that in the THz region the illumination conditions are non-homogeneous with an exponential decrease of illuminating signal strength towards the edges of the field-of-view. Due to the large coherence length in this region of the electromagnetic spectrum, one can not simply make the illuminating beam larger and use only the center portion for imaging. This will cause interference effects that appear as dominating artifacts in images (see Figure 4), which limit image quality and spatial resolution. So, without calibration one either has to live with a limited field-of-view or one has to accept image artifacts and limited spatial resolution. Therefore, calibration for non-homogeneous illumination is essential for a practical THz single-pixel imaging system. We will discuss this issue and possible solutions in the next section. Note that for a practical imaging system the calibration of the illuminating beam is very important and also leads to challenges related to the used masks. This is also exemplified in Figure 5, which shows a $0.35\text{\,}\mathrm{THz}$ single-pixel camera measurement of a non-metallic Siemenssstar test target reconstructed using a convolution approach with a non-negative least squares approach, i.e. we use convolutional masks leading to 2D circulant matrices.

In a modulated illumination setting of a THz single-pixel camera the quality and fidelity of the masks/illumination patterns determine the achievable spatial resolution, the signal-to-noise ratio, the achievable undersampling ratio and even more. The physical process of implementing the masks is, therefore, very important and potentially introduces deviations into the masks already on the physical level. These potential deviations can be simple blurring effects, the introduction of an offset or the reduction in the so-called modulation depth. As introduced, for binary masks the modulation depth is essentially the difference between ones and zeros. So, the ideal value is of course, 1 or 100% but in a practical system the modulation depth can be several 10 percent below the ideal value. Depending on the chosen reconstruction approach this will severely influence the image fidelity of reconstructed images (see Figure 5 for an example). The example in Figure 5 shows what happens when the modulation depth is only 40% and the masks are not optimized in an undersampling modality. Due to the fact that the spatial resolution and the signal-to-noise ratio in the image are severely limited the image appears blurry and noisy.

All the aforementioned issues have to be considered on the software level in order to harness the potential power of a THz single-pixel camera and, therefore, often robustness is an important consideration when choosing the reconstruction method. As mentioned the reconstruction approach used in this example was based on the idea of convolutional mask, which offers strong simplifications in the design of masks and the memory consumptions, moreover the reconstruction algorithms can be made more efficient. Hence, such an approach has a lot of potential for THz single-pixel cameras, but there is still a lot of effort necessary in order to optimize the imaging process.

5 Calibration Problems

In the previous section we have seen the difficulties to calibrate the illumination beam directly, hence it seems necessary to perform a self-calibration approach during the reconstruction. We hence look for a multiplicative decomposition of the image $\mathbf{p}$ into a smooth light field $\mathbf{d}$ and a normalized target structure $\mathbf{x}$, i.e.

[TABLE]

where $\odot$ denotes pointwise multiplication. A simple approach would be to first reconstruct the image and then use standard decomposition ideas. However, this approach seems suboptimal for undersampled measurements, since better a-priori knowledge on light and target is available than for their composition.

5.1 Self-Calibration and Bilinear Inverse Problems

A very recent trend in compressed sensing is to include the calibration task into the recovery step as well. This approach has been termed also as self-calibration. In the case of unknown illumination this yields a bilinear inverse problem. In particular, for sparsity promoting regularizers this falls in the category of biconvex compressed sensing [LS15]. In this work the lifting approach has been adopted to transform the task to an convex formulation, see here also the seminal work on Phaselift [CSV13]. Unfortunately, this approach does scale for imaging problems.

Here we are confronted with the generic problem:

[TABLE]

with an appropriate regularization functional $r$ for both structures, the illumination $\mathbf{d}$ and the target $\mathbf{x}$. This problem is a particular case of a bilinear inverse problem and linked to compressive blind deconvolution. Indeed, let us formalize this by using the 2D fast Fourier transform (FFT):

[TABLE]

where $\mathbf{F}\colon\mathbb{R}^{n}\rightarrow\mathbb{C}^{n}$ defines the Fourier transform operator and $\mathbf{l},\mathbf{r}\in\mathbb{C}^{n}$. and therefore the observation $\mathbf{y}=\frac{1}{\sqrt{n}}\mathbf{A}\mathbf{F}^{-1}(\mathbf{l}\ast\mathbf{r})$ is a compressed circular 2D convolution of $\mathbf{l}$ and $\mathbf{r}$. Obviously, without further assumptions this type of inverse problem cannot be solved uniquely.

Let us discuss the case $\mathbf{A}=\text{Id}$ that usually corresponds to an image scanning approach (and not the single pixel setup). Here, the measurement image is $\mathbf{p}$ and the goal is to factorize it into light $\mathbf{d}$ and a normalized target $\mathbf{x}$ using the program:

[TABLE]

Note that this formulation is a non-convex problem since it is not jointly convex in $(\mathbf{d},\mathbf{x})$. To make this problem more well-posed from a compressed viewpoint further assumptions on $\mathbf{d}$ and $\mathbf{x}$ are necessary. In the case of random subspace assumptions for either $\mathbf{d}$ or $\mathbf{x}$ and some additional incoherence conditions, recovery guarantees for the convexified (lifted) formulation of blind deconvolution have been established in [ARR14]. This framework even extends to the more involved problem of blindly demixing convolutional mixtures, called as also blind demixing and deconvolution, and recently almost-optimal sampling rates could established here under similar assumptions [JKS18]. However, the analysis of the original non-convex problem (without lifting) is often difficult, requires a priori assumptions and sufficient randomization and the performance of iterative decent algorithms often depends on a good initialization. First results in this direction appeared here for example in [LLSW18].

A further way to convexify a related problem is based on a formulation for strictly positive gains

[TABLE]

In this problem the image $\mathbf{x}$ is unknown and each measurement has an unknown strictly-positive gain $\mathbf{g}>\epsilon$. This problem occurs exactly in our setting when $\mathbf{A}=\text{Id}$ and then $\mathbf{g}=\mathbf{d}$. Thus, under the additional prerequiste of $\mathbf{d}\geq\epsilon>0$ one could write a constraint $\mathbf{p}=\mathbf{d}\odot\mathbf{x}$ also as $\mathbf{d}^{-1}\odot\mathbf{p}=\mathbf{x}$ which is jointly convex in $(\mathbf{d}^{-1},\mathbf{x})$ [GCD12] [LS18]. This approach is also interesting if $\mathbf{x}=\mathbf{\Psi}\mathbf{\boldsymbol{\alpha}}$ itself is a compressed representation. Unfortunately, this approach does not apply to the single pixel imaging setup due to matrix $\mathbf{A}$, which cannot be interchanged with the division by $\mathbf{d}$.

5.2 Single Pixel and Retinex Denoising

A related classical problem in imaging is the so-called Retinex approach. The corresponding problem was first investigated by Land (cf. [Lan64] [LM71]) in the context of human visual studies and provide a first entry in the developing Retinex theory. The human eye is able to compensate lack of illumination, hence it is able to filter for structural information in the field of view. From a more general perspective this translates to the question of how to separate an image into a structural and an illumination part. Here we focus on the approach of Kimmel et al. [KES*+*03] and further developments of this approach. By defining $\mathbf{s}\mathrel{\mathop{:}\!\!=}\log(\mathbf{p})$, $\mathbf{r}\mathrel{\mathop{:}\!\!=}\log(\mathbf{x})$ and $\mathbf{l}\mathrel{\mathop{:}\!\!=}\log(\mathbf{d})$ we move the image into the logarithmic domain. Again in the case $\mathbf{A}=\text{Id}$ this allows for the usage of a convex variational model as proposed in [KES*+*03]. The basic assumptions are as follows:

• •

spatial smoothness of illumination,

• •

$\mathbf{l}$* is greater than $\mathbf{s}$*: Since $\mathbf{x}\in[0,1]^{n}$ and $\mathbf{d}\geq\mathbf{x}\odot\mathbf{d}$ it is $\mathbf{l}\geq\mathbf{s}$ due to the monotone nature of the logarithmic map,

• •

non trivial illumination: resemblance to the original image, hence $l-s$ small,

• •

soft smoothing on structure reconstruction and

• •

smooth boundary condition.

These result in the variation approach proposed by Kimmel et al.:

[TABLE]

However, this applies a smoothing on $l$ and $l-s$, which is basically $r$. Then this boils down to the data fitting term and two smoothing regularization operators that are balanced using the respective parameters $\alpha$ and $\beta$. In a later paper by Ng and Wang [NW11], the usage of TV regularization for the structural image was introduced. The proposed model is

[TABLE]

This equation makes sense from a Retinex point of view. One has given an image $\mathbf{p}$ or $\log(\mathbf{p})=\mathbf{s}$ respectively, and wants to separate reflection $\log(\mathbf{x})=\mathbf{r}$ and $\log(\mathbf{d})=\mathbf{l}$.

A comparison of both Retinex models are shown in Figure 7, with an interesting parameter dependent behavior. The regularization and weighting parameters $\alpha$ and $\beta$ have been determined using a parameter test, shown in Figure 6. We see that choosing parameters optimizing the SSIM measure for the target respectively its logarithm $\mathbf{r}$ ($\alpha=10^{0}$, $\beta=10^{2}$ in the Kimmel model respectively $\alpha=10^{1}$, $\beta=10^{2}$ for the Ng and Wang model) yields a strange result with respect to the illumination, which still has quite some of the target structure included, in addition to the full beam. This can be balanced by different parameter choices ($\alpha=10^{0}$, $\beta=10^{3}$ in the Kimmel model respectively $\alpha=10^{1}$, $\beta=10^{3}$ for the Ng and Wang model) displayed in the lower line respectively, which eliminates most structure from the reconstructed illumination, but also leaves some beam effects in the target.

However in the framework of single pixel camera approaches, one does not have the fully recovered image at hand. Instead we only have measured data, hence the data fidelity term in (15) is replaced by a reconstruction based data fidelity term based on the forward operator $A$ and measured data $y$:

[TABLE]

A simple approach is to apply a two-step idea: In the first step one can compute standard reconstruction of $\mathbf{p}$ e.g. by (6) and subsequently apply the Retinex model with $\mathbf{s}=\log\mathbf{p}$. Since in this case the first reconstruction step does not use prior knowledge about the structure of illumination, it is to be expected that the results are worse compared to a joint reconstruction when the number of measurements decreases.

In order to construct a computational approach for the calibration problem we still formulate the problem in the logarithmic variables

[TABLE]

and derive a forward-backward splitting algorithm in these variables: The forward step is simply given by

[TABLE]

and the backward step then computes $\mathbf{r}^{k+1}$ as a minimizer of

[TABLE]

and $\mathbf{l}^{k+1}$ as a minimizer of

[TABLE]

Note that by adding the two equations in the forward step we can directly formulate it as an update for $\mathbf{s}$ as

[TABLE]

This induces a simple idea for the two-step approach: we can iterate $\mathbf{s}$ with the above scheme and directly apply the above Retinex model to data $\mathbf{s}^{k+1/2}$. From the resulting minimizers $\mathbf{r}^{k+1}$ and $\mathbf{l}^{k+1}$ we can compute $\mathbf{s}^{k+1}=\mathbf{r}^{k+1}+\mathbf{l}^{k+1}$. The resulting reconstructions are shown in Figure 8, which are clearly suboptimal since too much of the random structure from the matrix $\mathbf{A}$ is propagated into $\mathbf{s}^{k+1/2}$, which is not reduced enough in the Retinex step. The results of the forward-backward splitting approach are shown in Figure 9, which are clearly improving the separation of illumination beam and target and yield robustness with respect to undersampling, although still not providing perfect smoothing of the images. This might be expected from improved regularization models for the decomposition to be investigated in the future.

6 Phase Retrieval in Single Pixel Cameras

In the following we discuss another aspect of reconstruction in single pixel cameras, namely phase reconstruction problems caused by diffraction effects. This problem is a particular instance of the difficult phase retrieval problem. Compared to blind deconvolution, as a (non-convex) bilinear inverse problem, phase retrieval is the corresponding quadratic (non-convex) case and therefore most analytical results here are based on lifting the problem to a convex formulation [CSV13]. For some overview, further references and due to limited space we refer here exemplary to the overview article [JEH16]. Interestingly also, first works already appeared where phase retrieval and blind deconvolution are combined using lifting [AAH18]. However, already for imaging problems at moderate resolution these approaches usually not scale and require a priori random subspace assumptions which are difficult to fulfill in practice. In the following we will discuss how to setup the phase reconstruction problem in the diffraction case.

The propagation of light waves can - after some approximation - be represented by the Fresnel diffraction integral, we refer to [Nic18] and references cited therein for a detailed treatment of the discretization problem. In the single pixel setup (compare Figure 1) we consider that the diffraction has to be taken into account in two ways, namely between the object and mask plane and the mask and detector plane. This means that the measurement matrix is further changed by the introduction of the masks. In addition to this only the magnitude of the combined complex signals is obtained. Let ${\bf D}_{om}$ be a matrix discretizing the diffraction integral from object to mask plane and ${\bf D}_{md}$ be a matrix discretizing the diffraction from mask to detector plane. Then the complex signal arriving at the detector plane is

[TABLE]

where $\mathbf{p}$ is the complex image, $\mathbf{A}_{i}$ denotes the $i$-th row of $\mathbf{A}$, and ${\bf 1}$ a vector filled with all entries equal to one. Note that in the absence of diffraction, i.e. ${\bf D}_{md}={\bf D}_{om}={\bf Id}$, this reduces to the standard single pixel camera approach discussed above. The measured intensity is then the absolute value of ${\bf z}_{i}$, or rather its square, i.e.

[TABLE]

Defining a matrix ${\bf B}$ with rows

[TABLE]

the phase reconstruction problem can be written in compact notation as

[TABLE]

It is apparent that (18) is a nonlinear problem compared to the linear phase- and diffractionless problem, where both the matrix $\mathbf{A}$ and the image $\mathbf{p}$ can be modeled to be real non-negative (hence the absolute value does not lead to a loss of information). For this reason the iterative solution of (18) respectively least-squares versions thereof is a problem of central importance, which we will discuss in the next section.

6.1 Algorithms for Phase Retrieval

Phase retrieval is a classical problem in signal processing, which has been revived by compressed sensing approaches in the last decade, cf. [Fie82] [BCL02] [Mar07] [JEH15] for an overview of classical and recent methods. A famous early algorithm is the Gerchberg-Saxton (GS) algorithm (cf. [Ger72]). However a version of the GS algorithm has been developed by Fienup (cf. [Fie82]) that works with amplitude measurements only and is using a multiplicative update scheme of the form

[TABLE]

where ${\bf B}^{+}$ defines the commonly used pseudo inverse of ${\bf B}$ (cf. [EHN00]). In the case of additional constraints on $\mathbf{p}$, those are applied to modify $\mathbf{y}_{k+1}$ in an additional projection step. Note that the original Gerchberg-Saxton algorithm is formulated for Fourier measurements only, where ${\bf B}$ is invertible by its adjoint. The algorithm can be used in particular to compute real images, where

[TABLE]

Another approach that can be found at several instances in literature (cf. [BKE97] [SBE14] [GX17]) is based on the application of (variants of) Gauss-Newton methods to the least squares problem of minimizing

[TABLE]

This amounts to linearizing the residuals around the last iterate and to obtain $\mathbf{p}_{k+1}$ as the minimizer of

[TABLE]

Several variants are used to stabilize the Gauss-Newton iteration and to account for the ill-conditioning of the Jacobian matrix

[TABLE]

A popular one, which we will also use in our numerical tests, is the Levenberg-Marquardt method, which computes $\mathbf{p}_{k+1}$ as the minimizer of

[TABLE]

with $\alpha_{k}$ a decreasing sequence of positive parameters.

A recently popular approach to solve the nonconvex least-squares problem is the Wirtinger flow (cf. [CLS15]). Its main ingredient is just gradient descent on the least squares functional $L$, i.e.,

[TABLE]

For the choice of the step size an increasing strategy like $\mu_{k}\sim 1-e^{-k/k_{0}}$ is proposed. Obviously the gradient descent as well as the Levenberg-Marquardt method above rely on appropriate initializations in order to converge to a global minimum. For this sake a spectral initialization has been proposed, which chooses $p_{0}$ as the first eigenvector of the positive semidefinite matrix

[TABLE]

corresponding to the data sensitivity in the least squares functional. For practical purposes the first eigenvector can be approximated well with the power method.

A very recent approach is the truncated amplitude flow (cf. [WGE18]), whose building block is gradient descent on the alternative least squares functional

[TABLE]

which is however not differentiable for ${\bf B}_{i}^{T}\mathbf{p}=0$. In the case the derivative exists we find

[TABLE]

The truncated amplitude flow now only selects a part of the gradient in order to avoid the use of small $|{\bf B}_{i}^{T}\mathbf{p}|$ (compared to $\sqrt{\mathbf{y}_{i}}$), i.e.

[TABLE]

with index set

[TABLE]

with some $\gamma>0$. For the truncated amplitude flow another initialization strategy has been proposed in [WGE18], which tries to minimize the angles of $\mathbf{p}_{0}$ to the measurement vectors ${\bf B}_{i}$ for a subset of indices $I_{0}$. This is equivalent to computing the eigenvector for the smallest eigenvalue of the matrix

[TABLE]

Since computing such an eigenvector is a problem of potentially high effort, it is approximated by computing the largest eigenvalue of a matrix built of the remaining rows of ${\bf B}$ (cf. [WGE18] for further details).

In order to introduce some regularization into any of the flows above, we can employ a forward-backward splitting approach. E.g. we can produce the Wirtinger flow to produce the forward estimate and use an additional backward proximal step

[TABLE]

with the regularization functional $R$ and a parameter $\lambda_{k+1}$ chosen appropriately in dependence of $\mu_{k+1}$. The proximal map $\mathbf{p}=\text{prox}_{\lambda_{k+1}R}(\mathbf{p}_{k+1/2})$ is given as the unique minimizer of

[TABLE]

Finally we mention that in addition to the nonconvex minimization approaches there has been a celebrated development towards convexifying the problem, the so-called PhaseLift approach, which can be shown to yield an exact convex relaxation under appropriate conditions (cf. [CSV13]). The key idea is to embed the problem into the space of $n\times n$ matrices and to find $\mathbf{p}\mathbf{p}^{T}$ as a low rank solution of a semidefinite optimization problem. The squared absolute value is rewritten as

[TABLE]

with the positive semidefinite rank-one matrix ${\bf P}=\mathbf{p}\mathbf{p}^{*}$. The system of quadratic equations $|{\bf B}_{i}\mathbf{p}|^{2}=y_{i}$ then is rewritten as the problem of finding the semidefinite matrix of lowest rank that solves the linear matrix equations ${\bf B}_{i}^{*}{\bf P}{\bf B}_{i}=y_{i}$. Under appropriate condition this problem can be exactly relaxed to minimizing the nuclear norm of ${\bf P}$, a convex functional, subject to linear equation and positive semidefiniteness. The complexity of solving the semidefinite problem in dimension $n^{2}$ makes this aprroach prohibitive for applications with large $n$ however, hence we will not use it in our computational study in the next section.

6.2 Results

In the following we present some results obtained from the algorithms discussed above when applied to our single pixel setup and compare their performance. We start by using the algorithms from the previous section for reconstructing a real-valued image showing the figure $\pi$ ($n=900$) without noise in the measurements. For this sake we consider different types of undersampling and three different initializations: a random initial value (rand), the spectral initialization proposed originally for the Wirtinger flow (spec) and the orthogonality promoting initialization originally proposed for the truncated amplitude flow (amp). The results for the four different methods are shown in Figures 10 and 11. The results clearly demonstrate the dependence of solutions on the initial value, in particular we see the improvement of the two new initialization strategies with respect to the random one (except in the case of the algorithm by Fienup, which however yields inferior results in all cases). We observe that the Levenberg-Marquardt method performs particularly well for larger sampling rates and yields an almost perfect reconstruction, but also produces suitable reconstructions for stronger undersampling. A more quantitative evaluation is given in Figure 12, which provides plots of the relative mean-square error vs. the sampling rate, surprisingly the Levenberg-Marquardt method outperforms the other schemes for most rates.

The positive effect of regularization and its necessity for very strong undersampling is demonstrated in Figure 13, where we display the reconstruction results obtained with the forward-backward splitting strategy and the total variation as regularization functional. We see that both approaches yield similar results and in particular allow to proceed to very strong undersampling with good quality results. The effect of noise in the data, here for a fixed sampling ratio $\frac{m}{n}=0.7$ is demonstrated in Figure 14 for different signal-to-noise ratios. Again, not surprisingly, the regularized version of the flows yields stable reconstructions of good quality even for data of lower quality.

Let us finally turn our attention to the reconstruction of the phase in a complex image. The ground truth for the amplitude and phase of the complex signal are shown in Figure 15. For brevity we only provide visualizations of reconstructions for the truncated amplitude flow and the Levenberg-Marquardt method in Figure 16, which clearly indicates that the truncated amplitude flow outperforms the Levenberg-Marquardt method, which is the case for all conducted experiments. Again, a more quantitative evaluation is given in Figure 17, which provides a plot of the relative mean-square error (made invariant to a global phase that is not identifiable) vs. the sampling rate. We see that the Wirtinger flow can obtain the same reconstruction quality for very low sampling rates, but is outperformed by the truncated amplitude flow at higher sampling rates.

7 Conclusion

We have seen that the compressed sensing approach based on single-pixel imaging has great potential to decrease measurement time and effort in THz imaging, but the application in practice depends on several further challenges to be mastered. First of all appropriate regularization models and numerical algorithms are needed for the image reconstruction problem in order to obtain higher image resolution at reasonable computational times. Moreover, in several situations it is crucial to consider auto-calibration issues in particular related to the fact that the illuminating beam is difficult to be characterized, and in some cases also diffraction becomes relevant, which effectively yields a phase retrieval problem. Both effects change the image reconstruction from a problem with a rather simple and incoherent linear forward model to a nonlinear problem with a more coherent forward model, which raises novel computational and also theoretical issues, since the assumptions of the existing compressed sensing theory are not met.

Besides the above mentioned issues several further aspects of practical imaging are foreseen to become relevant for THz single pixel imaging, examples being motion correction for imaging in the field or of moving targets and the reconstruction of multi- or hyper-spectral images, which is a natural motivation in the THz regime. In the latter case it will become a central question how to combine as few masks as possible in different spectral locations.

Acknowledgement

This work has been supported by the German research foundation (DFG) through SPP 1798, Compressed Sensing in Information Processing, projects BU 2327/4-1 and JU 2795/3. MB acknowledges further support by ERC via Grant EU FP7 ERC Consolidator Grant 615216 LifeInverse.