Reconstruction Methods in THz Single-pixel Imaging

03/21/2019 ∙ by Martin Burger, et al. ∙ 0

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.



There are no comments yet.


page 8

page 10

page 14

page 15

page 20

page 21

page 23

page 25

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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. [D06]), 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. [DDT08] [CCT08] [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. [CCT08] [WSM14] [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

Figure 1: Sketch of the setup of imaging planes in THz single pixel imaging (from [Nic18]).

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 , i.e. where . The simplest model of the image formation process is to have each measurement , , as the sum of subpixel values for those corresponding to the open subpixels in the -th mask. This means we can write


with a matrix whose entries are only zeros and ones. Again the non-zero entries in each row of correspond to the open subpixels of the respective mask (compare Figure 1 for a schematic overview of such a setup).

Choosing deterministic masks appropriately one could guarantee that is invertible and simply solve the linear reconstruction problem. However, in practice one would like to reconstruct a reasonable image from a number of measurements significantly smaller than . 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

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 images are compressible and can be (approximately) sparsely represented as with and a particular basis , e.g. wavelets, or overcomplete systems with such as shearlets (cf. [GL07]). By this we mean that


for small , is considerably smaller then the ambient dimension . Ideally, one would thus solve the problem of minimizing subject to the linear constraint , 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 the convex relaxation


recovers exactly (in the noiseless setting) the unknown yielding the correct image . A common additional assumption is that the image itself has non-negative pixel intensities, such that problem (2) is extended to:


If and if the row span of

intersects the positive orthant, meaning that there exists a vector

such that , the measurement matrix itself already assures that the -norm of a feasible in (3) is equal (or close in the noisy case) to where is the unknown image to recover. To see this, let us assume exemplary that we find a vector such that is the all-one vector and . Then

and hence, it is enough to minimize the residual over and replace (3) by a simple non-negative least squares (NNLS) problem:


Indeed, that under these assumptions -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 .

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


where is an appropriate regularization parameter.

2.2 Total Variation Regularization and Related Methods

While simple -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


or the penalty version


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 -norm of is used in the regularization (namely ), 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 (the rows corresponding to partial derivatives) and computes 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])


where is a parameter to be chosen appropriately. Again in practice suitable isotropic variants of the -norm are used, and in most approaches only the symmetrized gradient of the vector field 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

with being an appropriately chosen seminorm and a suitable matrix , e.g. the gradient or in the case of using a basis, or

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 as the minimizer of

with subgradient . In this way in each iteration a suitable distance to the last iterate is penalized instead of the original functional , which is effectively a distance to zero. The parameter is chosen much larger than the optimal in the above problem, roughly speaking when carrying out iterations we have Observing from the optimality condition that with , the Bregman iteration can be interpreted equivalently as the augmented Lagrangian method for the constrained problem of minimizing subject to , i.e.


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 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. has no particular structure. Thus, the typical source condition for some 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 -norm, whose -proximal map

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

with convex functionals and and a linear operator . The primal dual approach reformulates the minimization as a saddle point problem

with being the convex conjugate of . Then one iterates primal minimization with respect to and dual maximization with respect to 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.

Figure 2: Grayscale images of synthetic data, synthetic illumination beam with values between 0 (black) and 1 (white); top left: Gaussian beam; top center: x phantom ground truth; bottom center: x phantom with applied beam; top right: gradient phantom ground truth; bottom right: gradient phantom with applied beam.

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, - 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.

Figure 3: Reconstruction the x and gradient phantom for multiple sampling rates (rows) and reconstruction frameworks (columns).

4 Challenges in Practical THz Single-Pixel Imaging

Figure 4: Mechanically scanned image of a metal Siemensstar test target with a diameter of . The image acquisition took several hours () depending on the number of steps but still scanning artifacts are very prominent. The inset shows a photo of the metal Siemensstar target.

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 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.

Figure 5: Megapixel image using a convolution approach (i.e. a binary circulant matrix). Labels show the pixel numbers in both directions, about 10 pixel correspond to mm. The spatial resolution in the image is still limited by sub-optimal illumination patterns but the image shows that the resolution is reasonable using the convolution approach, while being very fast (more than one FPS can be achieved). The inset shows again a photo of the imaged object.

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 into a smooth light field and a normalized target structure , i.e.


where 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:


with an appropriate regularization functional for both structures, the illumination and the target . 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):


where defines the Fourier transform operator and . and therefore the observation is a compressed circular 2D convolution of and . Obviously, without further assumptions this type of inverse problem cannot be solved uniquely.

Let us discuss the case that usually corresponds to an image scanning approach (and not the single pixel setup). Here, the measurement image is and the goal is to factorize it into light and a normalized target using the program:


Note that this formulation is a non-convex problem since it is not jointly convex in . To make this problem more well-posed from a compressed viewpoint further assumptions on and are necessary. In the case of random subspace assumptions for either or 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


In this problem the image is unknown and each measurement has an unknown strictly-positive gain . This problem occurs exactly in our setting when and then . Thus, under the additional prerequiste of one could write a constraint also as which is jointly convex in [GCD12] [LS18]. This approach is also interesting if itself is a compressed representation. Unfortunately, this approach does not apply to the single pixel imaging setup due to matrix , which cannot be interchanged with the division by .

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. [KES03] and further developments of this approach. By defining , and we move the image into the logarithmic domain. Again in the case this allows for the usage of a convex variational model as proposed in [KES03]. The basic assumptions are as follows:

  • spatial smoothness of illumination,

  • is greater than : Since and it is due to the monotone nature of the logarithmic map,

  • non trivial illumination: resemblance to the original image, hence small,

  • soft smoothing on structure reconstruction and

  • smooth boundary condition.

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


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


This equation makes sense from a Retinex point of view. One has given an image or respectively, and wants to separate reflection and .

A comparison of both Retinex models are shown in Figure 8, with an interesting parameter dependent behavior. The regularization and weighting parameters and 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 (, in the Kimmel model respectively , 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 (, in the Kimmel model respectively , 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.

L2 based (see (14)) TV based (see (15))
Figure 6: Overview of PSNR value and SSIM index for systematic parameter test for the synthetic x phantom dataset. Map of regularization parameters in respect to given pair of regularization parameters and : in axes with values (top), , …, , (bottom) and in axis with values (left), , …, , (right).
Figure 7: Comparison of the L2 based Retinex model by Kimmel and the total variation based Retinex model by Ng and Wang for both datasets.
Figure 8: Reconstruction results for an alternating reconstruction and Retinex scheme.

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 and measured data :


A simple approach is to apply a two-step idea: In the first step one can compute standard reconstruction of e.g. by (6) and subsequently apply the Retinex model with . 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.

Figure 9: Reconstruction results for Forward-Backward Splitting including the L1 based model proposed by Ng and Wang.

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


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

and the backward step then computes as a minimizer of

and as a minimizer of

Note that by adding the two equations in the forward step we can directly formulate it as an update for as

This induces a simple idea for the two-step approach: we can iterate with the above scheme and directly apply the above Retinex model to data . From the resulting minimizers and we can compute . The resulting reconstructions are shown in Figure 8, which are clearly suboptimal since too much of the random structure from the matrix is propagated into , 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 be a matrix discretizing the diffraction integral from object to mask plane and be a matrix discretizing the diffraction from mask to detector plane. Then the complex signal arriving at the detector plane is

where is the complex image, denotes the -th row of , and a vector filled with all entries equal to one. Note that in the absence of diffraction, i.e. , this reduces to the standard single pixel camera approach discussed above. The measured intensity is then the absolute value of , or rather its square, i.e.

Defining a matrix with rows

the phase reconstruction problem can be written in compact notation as


It is apparent that (18) is a nonlinear problem compared to the linear phase- and diffractionless problem, where both the matrix and the image 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

where defines the commonly used pseudo inverse of (cf. [EHN00]). In the case of additional constraints on , those are applied to modify in an additional projection step. Note that the original Gerchberg-Saxton algorithm is formulated for Fourier measurements only, where is invertible by its adjoint. The algorithm can be used in particular to compute real images, where

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


This amounts to linearizing the residuals around the last iterate and to obtain as the minimizer of


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

A popular one, which we will also use in our numerical tests, is the Levenberg-Marquardt method, which computes as the minimizer of


with 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 , i.e.,


For the choice of the step size an increasing strategy like 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

as the first eigenvector of the positive semidefinite matrix

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


which is however not differentiable for . In the case the derivative exists we find

The truncated amplitude flow now only selects a part of the gradient in order to avoid the use of small (compared to ), i.e.


with index set

with some . For the truncated amplitude flow another initialization strategy has been proposed in [WGE18], which tries to minimize the angles of to the measurement vectors for a subset of indices

. This is equivalent to computing the eigenvector for the smallest eigenvalue of the matrix

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 (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


with the regularization functional and a parameter chosen appropriately in dependence of . The proximal map is given as the unique minimizer of

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 matrices and to find as a low rank solution of a semidefinite optimization problem. The squared absolute value is rewritten as

with the positive semidefinite rank-one matrix . The system of quadratic equations then is rewritten as the problem of finding the semidefinite matrix of lowest rank that solves the linear matrix equations . Under appropriate condition this problem can be exactly relaxed to minimizing the nuclear norm of , a convex functional, subject to linear equation and positive semidefiniteness. The complexity of solving the semidefinite problem in dimension makes this aprroach prohibitive for applications with large however, hence we will not use it in our computational study in the next section.

6.2 Results

Figure 10: Reconstructions of -image with the Fienup-variant of the Gerchberg-Saxton algorithm and the Wirtinger flow (from [Nic18]).
Figure 11: Reconstructions of -image with the truncated amplitude flow and the Levenberg-Marquardt method (from [Nic18]).
Figure 12: Plot of relative mean square error to ground truth vs. sampling rate for reconstructions of -image (from [Nic18]).

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 () 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 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.

Figure 13: Reconstructions of -image with the regularized version of the truncated amplitude flow and the Wirtinger flow (from [Nic18]).
Figure 14: Reconstructions of -image from noisy data with different signal-to-noise ratios (SNR) (from [Nic18]).

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.

Figure 15: Ground truth for amplitude (left) and phase (right, from [Nic18]).
Figure 16: Reconstructions of amplitude (left) and phase (right, from [Nic18]) with truncated amplitude flow and Levenberg-Marquardt method with different initializations and sampling rates.
Figure 17: Relative mean-square error in the reconstructions of the complex image (from [Nic18]).

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.


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.


  • [AAH18] Ali Ahmed, Alireza Aghasi, and Paul Hand. Blind deconvolutional phase retrieval via convex programming. CoRR, 2018.
  • [AFJH17] S Augustin, S Frohmann, P Jung, and H-W Hübers. An optically controllable 0.35 thz single-pixel camera for millimeter resolution imaging. In 2017 42nd International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz), pages 1–2. IEEE, 2017.
  • [AHJH15] Sven Augustin, J Hieronymus, P Jung, and H-W Hübers. Compressed sensing in a fully non-mechanical 350 ghz imaging setting. Journal of Infrared, Millimeter, and Terahertz Waves, 36(5):496–512, 2015.
  • [ARR14] Ali Ahmed, Justin Romberg, and Benjamin Recht. Blind deconvolution using convex programming. IEEE Trans. Inf. Theory, 60(3):1711–1732, 2014.
  • [Bar07] R. G. Baraniuk. Compressive sensing. IEEE Signal Processing Magazine, 24(4):118–121, July 2007.
  • [BB18] Martin Benning and Martin Burger. Modern regularization methods for inverse problems. Acta Numerica, 27:1–111, 2018.
  • [BCL02] Heinz H Bauschke, Patrick L Combettes, and D Russell Luke. Phase retrieval, error reduction algorithm, and fienup variants: a view from convex optimization. JOSA A, 19(7):1334–1345, 2002.
  • [BDS18] Martin Burger, Hendrik Dirks, and Carola-Bibiane Schonlieb. A variational model for joint motion estimation and image reconstruction. SIAM Journal on Imaging Sciences, 11(1):94–128, 2018.
  • [BEZ08] Alfred M Bruckstein, Michael Elad, and Michael Zibulevsky. On the uniqueness of non-negative sparse & redundant representations. In IEEE International Conference on Acoustics, Speech and Signal Processing, pages 5145–5148. IEEE, 2008.
  • [BKE97] Barbara Blaschke-Kaltenbacher and Heinz W Engl. Regularization methods for nonlinear ill-posed problems with applications to phase reconstruction. In Inverse Problems in Medical Imaging and Nondestructive Testing, pages 17–35. Springer, 1997.
  • [BKP10] Kristian Bredies, Karl Kunisch, and Thomas Pock. Total generalized variation. SIAM Journal on Imaging Sciences, 3(3):492–526, 2010.
  • [BO13] Martin Burger and Stanley Osher. A guide to the tv zoo. In Level set and PDE based reconstruction methods in imaging, pages 1–70. Springer, 2013.
  • [BSS16] Martin Burger, Alexander Sawatzky, and Gabriele Steidl. First order algorithms in variational image processing. In Splitting Methods in Communication, Imaging, Science, and Engineering, pages 345–407. Springer, 2016.
  • [CCN14] Vicent Caselles, Antonin Chambolle, and Matteo Novaga. Total variation in imaging. Handbook of Mathematical Methods in Imaging, pages 1–39, 2014.
  • [CCT08] Wai Lam Chan, Kriti Charan, Dharmpal Takhar, Kevin F Kelly, Richard G Baraniuk, and Daniel M Mittleman. A single-pixel terahertz imaging system based on compressed sensing. Applied Physics Letters, 93(12):121105, 2008.
  • [CLS15] Emmanuel J Candes, Xiaodong Li, and Mahdi Soltanolkotabi. Phase retrieval via wirtinger flow: Theory and algorithms. IEEE Transactions on Information Theory, 61(4):1985–2007, 2015.
  • [CP11a] Emmanuel J Candes and Yaniv Plan. A probabilistic and ripless theory of compressed sensing. IEEE transactions on information theory, 57(11):7235–7254, 2011.
  • [CP11b] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of mathematical imaging and vision, 40(1):120–145, 2011.
  • [CRT06] E. J. Candes, J. Romberg, and T. Tao. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inform. Theory, 52(2):489–509, 2006.
  • [CSV13] Emmanuel J Candes, Thomas Strohmer, and Vladislav Voroninski. Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming. Communications on Pure and Applied Mathematics, 66(8):1241–1274, 2013.
  • [CX15] Jian-Feng Cai and Weiyu Xu. Guarantees of total variation minimization for signal recovery. Information and Inference, pages 328–353, may 2015.
  • [D06] David L Donoho et al. Compressed sensing. IEEE Transactions on information theory, 52(4):1289–1306, 2006.
  • [DDT08] Marco F Duarte, Mark A Davenport, Dharmpal Takhar, Jason N Laska, Ting Sun, Kevin F Kelly, and Richard G Baraniuk. Single-pixel imaging via compressive sampling. IEEE signal processing magazine, 25(2):83–91, 2008.
  • [Dir16] Hendrik Dirks. A flexible primal-dual toolbox. arXiv:1603.05835 [cs, math], 2016.
  • [Don06] D. L. Donoho. Compressed sensing. IEEE Trans. Inform. Theory, 52(4):1289–1306, 2006.
  • [EGP18] Matthew P Edgar, Graham M Gibson, and Miles J Padgett. Principles and prospects for single-pixel imaging. Nature Photonics, page 1, 2018.
  • [EHN00] Heinz Werner Engl, Martin Hanke, and Andreas Neubauer. Regularization of Inverse Problems, volume 375. Springer Netherlands, 1 edition, 2000.
  • [EK12] Yonina C Eldar and Gitta Kutyniok. Compressed sensing: theory and applications. Cambridge University Press, 2012.
  • [Fie82] James R Fienup. Phase retrieval algorithms: a comparison. Applied optics, 21(15):2758–2769, 1982.
  • [FR13] S. Foucart and H. Rauhut. A Mathematical Introduction to Compressive Sensing. Applied and Numerical Harmonic Analysis. Birkhäuser/Springer, New York, 2013.
  • [GCD12] Rémi Gribonval, Gilles Chardon, and Laurent Daudet. Blind calibration for compressed sensing by convex optimization. In 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 2713–2716. IEEE, 2012.
  • [Ger72] Ralph W Gerchberg. A practical algorithm for the determination of phase from image and diffraction plane pictures. Optik, 35:237–246, 1972.
  • [GL07] Kanghui Guo and Demetrio Labate. Optimally sparse multidimensional representation using shearlets. SIAM journal on mathematical analysis, 39(1):298–318, 2007.
  • [GX17] Bing Gao and Zhiqiang Xu. Phaseless recovery using the Gauss–Newton method. IEEE Transactions on Signal Processing, 65(22):5885–5896, 2017.
  • [JEH15] Kishore Jaganathan, Yonina C Eldar, and Babak Hassibi. Phase retrieval: An overview of recent developments. arXiv preprint arXiv:1510.07713, 2015.
  • [JEH16] Kishore Jaganathan, Yonina C. Eldar, and Babak Hassibi. Phase Retrieval: An Overview of Recent Developments. In A. Stern, editor, Optical Compressive Imaging, pages 263–287. CRC Press, 2016.
  • [JKM18] Peter Jung, Richard Kueng, and Dustin G. Mixon. Derandomizing compressed sensing with combinatorial design. CoRR, 2018.
  • [JKS18] Peter Jung, Felix Krahmer, and Dominik Stoeger. Blind Demixing and Deconvolution at Near-Optimal Rate. IEEE Transactions on Information Theory, 64(2):704–727, feb 2018.
  • [KES03] Ron Kimmel, Michael Elad, Doron Shaked, Renato Keshet, and Irwin Sobel. A variational framework for retinex.

    International Journal of Computer Vision

    , 52(1):7–23, 2003.
  • [KJ17] Richard Kueng and Peter Jung. Robust Nonnegative Sparse Recovery and the Nullspace Property of 0/1 Measurements. IEEE Transactions on Information Theory, 64(2):689–703, 2017.
  • [KKS17] Felix Krahmer, Christian Kruschel, and Michael Sandbichler. Total variation minimization in compressed sensing. In Compressed Sensing and its Applications, pages 333–358. Springer, 2017.
  • [KL12] Gitta Kutyniok and Demetrio Labate. Shearlets: Multiscale analysis for multivariate data. Springer Science & Business Media, 2012.
  • [Lan64] Edwin H Land. The retinex. American Scientist, 52(2):247–264, 1964.
  • [LLKR15] Weizhi Lu, Weiyu Li, Kidiyo Kpalma, and Joseph Ronsin. Compressed sensing performance of random Bernoulli matrices with high compression ratio. IEEE Signal Processing Letters, 22(8):1074–1078, 2015.
  • [LLSW18] Xiaodong Li, Shuyang Ling, Thomas Strohmer, and Ke Wei. Rapid, robust, and reliable blind deconvolution via nonconvex optimization. Applied and Computational Harmonic Analysis, pages 1–49, 2018.
  • [LM71] Edwin H. Land and John J. McCann. Lightness and Retinex Theory. Journal of the Optical Society of America, 61(1):1, 1971.
  • [LS15] Shuyang Ling and Thomas Strohmer. Self-calibration and biconvex compressive sensing. Inverse Problems, 31(11), 2015.
  • [LS18] Shuyang Ling and Thomas Strohmer. Self-calibration and bilinear inverse problems via linear least squares. SIAM Journal on Imaging Sciences, 11(1):252–292, 2018.
  • [Mar07] Stefano Marchesini. A unified evaluation of iterative projection algorithms for phase retrieval. Review of scientific instruments, 78(1):011301, 2007.
  • [Nic18] Lukas Nickel. Phase retrieval in single detector cameras. Master’s thesis, WWU Münster, 2018.
  • [NW11] Michael K. Ng and Wei Wang. A total variation model for retinex. SIAM Journal on Imaging Sciences, 4(1):345–365, 2011.
  • [NW13a] Deanna Needell and Rachel Ward. Near-optimal compressed sensing guarantees for total variation minimization. IEEE transactions on image processing, 22(10):3941–3949, 2013.
  • [NW13b] Deanna Needell and Rachel Ward. Stable image reconstruction using total variation minimization. SIAM Journal on Imaging Sciences, 6(2):1035–1058, 2013.
  • [Poo15] Clarice Poon. On the role of total variation in compressed sensing. SIAM Journal on Imaging Sciences, 8(1):682–720, 2015.
  • [SBE14] Yoav Shechtman, Amir Beck, and Yonina C Eldar. Gespar: Efficient phase retrieval of sparse signals. IEEE transactions on signal processing, 62(4):928–938, 2014.
  • [SH13] Martin Slawski and Matthias Hein. Non-negative least squares for high-dimensional linear models: Consistency and sparse recovery without regularization. Electron. J. Statist., 7:3004–3056, 2013.
  • [SSPB18] Robert J Stokoe, Patrick A Stockton, Ali Pezeshki, and Randy A Bartels. Theory and applications of structured light single pixel imaging. In Three-Dimensional and Multidimensional Microscopy: Image Acquisition and Processing XXV, volume 10499, page 104990E. International Society for Optics and Photonics, 2018.
  • [WGE18] Gang Wang, Georgios B Giannakis, and Yonina C Eldar. Solving systems of random quadratic equations via truncated amplitude flow. IEEE Transactions on Information Theory, 64(2):773–794, 2018.
  • [WSM14] Claire M Watts, David Shrekenhamer, John Montoya, Guy Lipworth, John Hunt, Timothy Sleasman, Sanjay Krishna, David R Smith, and Willie J Padilla. Terahertz compressive imaging with metamaterial spatial light modulators. Nature Photonics, 8(8):605, 2014.
  • [WXT11] Meng Wang, Weiyu Xu, and Ao Tang. A unique ”nonnegative” solution to an underdetermined system: From vectors to matrices. IEEE Trans. Inform. Theory, 59(3):1007–1016, 2011.
  • [ZBO11] Xiaoqun Zhang, Martin Burger, and Stanley Osher. A unified primal-dual algorithm framework based on Bregman iteration. Journal of Scientific Computing, 46(1):20–46, 2011.