First order optimization methods for terahertz digital holography

• Max. student(s): 1
• Advisors: Laurent Jacques and Murielle Kirkove, signal processing group, CSL, ULiège
• Teaching Assistant: Olivier Leblanc

Imaging context:

Terahertz (THz) imaging uses sources and detectors that use long-wave electromagnetic radiation. These waves allow entry to a wide variety of materials that are opaque to the human eye. The Centre Spatial de Liège (CSL, Liège Space center) is developing coherent diffraction imaging methods that require image reconstruction calculations. One of these methods is digital holography (DH), which allows, on the basis of a simple image acquisition, to reconstruct the interior of an object crossed by a coherent radiation (from a laser). Few years ago, an algorithm method has been developed at the CSL to reconstruct images from THz observations.

Description:

The direct model of digital holography (DH) is an instance of the phase retrieval problem [KZG21, Fienup82], that is, the recovery of a signal or an image from the magnitude of its Fourier transform. Mathematically, the purpose of DH is to collect a set of (noisy) observations in a measurement map $\boldsymbol y$ explained by the following forward model:

\begin{equation} \boldsymbol y = |\boldsymbol A_{d} \boldsymbol \psi + \boldsymbol \alpha \odot \boldsymbol r|^{2} + \boldsymbol n \tag{1} \end{equation}

where \(\boldsymbol \psi\) is an unknown complex “image” representing the object of interest, \(\boldsymbol A_d\) is a (known) linear transformation (associated with the light propagation at distance dd ; a simple matrix if one represents \(\boldsymbol \psi\) as a vector), \(\boldsymbol r\) is a unit reference beam (a controlled oscillating function) used to probe the content of \(\boldsymbol \psi\) , the quantity \(\boldsymbol \alpha\) is an unknown map that modulates (by the element-wise product \(\odot\)) the reference beam ( \(\boldsymbol \alpha\) accounts for non-uniformity in the experimental setup), and nn represents an additive measurement noise corrupting the model. Note that in equation (1) the square modulus \(|\cdot|^2\) is applied individually to each component of the vector \(\boldsymbol A_{d} \boldsymbol \psi + \boldsymbol \alpha \odot \boldsymbol r\).

The objective of digital holography is thus to deduce from the recorded measurements \(\boldsymbol y\) (the data) the unknown complex image \(\boldsymbol \psi\) and the (unknown) map \(\boldsymbol \alpha\) ; in words, we solve the inverse problem related to the forward relation in equation (1).

Classical approaches solve (1) from “direct” calculations in the Fourier domain (see the figure above); in a nutshell, the reference beam is set to an oscillating function (experimentally, thanks to an off-axis setup) that shifts the spectral content of the image, and allows to recover it despite the phase loss in the initial measurements \(\boldsymbol y\) . However, these approaches are not robust to noise, they do not allow measurements subsampling (as in compressive sensing theory [Baraniuk07]), and they are often corrupted by spurious artifacts for unknown, non-constant modulations \(\boldsymbol \alpha\) of the reference field \(\boldsymbol r\) .

An alternative is therefore to find two estimates \(\hat{\boldsymbol \psi}\)​ and \(\hat{\boldsymbol \alpha}\) by minimizing the following fidelity term, or cost function, which measures the distance between the actual measurements and those defined by two candidate fields \(\{\boldsymbol \psi’, \boldsymbol \alpha’\}\):

\begin{equation} D(\boldsymbol \psi’, \boldsymbol \alpha’) := \| \boldsymbol y - | \boldsymbol A_{d} \boldsymbol \psi’ + \boldsymbol \alpha’ \odot \boldsymbol r|^{2} \|^{2}. \tag{2} \end{equation}

Compared to direct Fourier methods, this approach can potentially remove reconstruction artifacts, be more tolerant to noise, under-sampling and difficult recording conditions.

Solving this minimization is challenging. First, it is non-convex relatively to the unknown map \(\boldsymbol \psi\) and \(\boldsymbol \alpha\) ; therefore, minimizing \(D\) can lead to many non-optimal minima. Second, the problem can even be ill-posed: it can be undetermined—if the number of unknowns (that is, the dimensions of \(\boldsymbol \psi\) and \(\boldsymbol \alpha\) ) is greater than the number of measurements (the dimension of \(\boldsymbol y\) )—, highly sensitive to noise.

Therefore, it was proposed in [KZG21] to regularize this minimization to stabilize the problem and select plausible solutions. They thus proposed to estimate \(\hat{\boldsymbol \psi}\)​ and \(\hat{\boldsymbol \alpha}\) by minimizing this new cost function:

\begin{equation} D(\boldsymbol \psi’, \boldsymbol \alpha’) + R_{\psi}(\boldsymbol \psi’) + R_{\alpha}(\boldsymbol \alpha’). \tag{3} \end{equation}

In words, two new functions, \(R_{\psi}\) and \(R_{\alpha}\) are added to the previous cost to ensure certain structures or smoothness to the fields \(\{\boldsymbol \psi’, \boldsymbol \alpha’\}\). In practice, the proposed regularization is done in the wavelet transform domain, and the two previous functions are set to the \(\ell_{1}\)-norm of wavelet coefficients of both \(\boldsymbol \psi’\) and \(\boldsymbol \alpha’\). This promotes the sparsity of these two estimates; they are described only by a few parameters in the wavelet domains, which helps to solve the ill-posedness of the former minimization.

In [KZG21], the regularized optimization above is solved with a dedicated algorithm using an alternating direction method of multipliers (ADMM) based framework. The figure on the top-right shows an example of reconstruction on simulated data.

Master project objectives:

While promising, the method proposed in [KZG21] shows limitations: it yields blocky artifacts and has convergence problems for too strong regularizations. Consequently, it delivers a non-optimal separation between the solutions \(\boldsymbol \psi\) and \(\boldsymbol \alpha\) .

The objective of this master project will thus to propose novel algorithmic strategies to solve the inverse problem (1).

In particular, it will be asked to the interested student to:

1. study, understand and implement the forward model (1) and the cost function (2) (in Python and/or Matlab);
2. implement first order optimization methods (such as proximal algorithms [PB14]) to minimize (3);
3. study how one can accumulate observations (with various reference fields \(\boldsymbol r\) ) to improve the quality of the reconstructed images;
4. validate the designed methods on numerical simulations as well as on actual data recorded by the CSL in Liège.

Supervision:

The master thesis will be carried out under the supervision of Pr. Laurent Jacques and Olivier Leblanc (ICTEAM institute at UCLouvain), and Murielle Kirkove at the signal processing group of CSL, in collaboration with Y. Zhao in CSL, Liège, for the experimental acquisition. Regular meetings (every week or every other week) will be organized during the year. The results of the master thesis are expected to develop the field of THz digital holography reconstruction. Depending on the achieved results, the work can lead to the submission of a paper in this scientific field (in a conference or a workshop).

References:

• [Fienup82] Fienup, J. R. (1982). Phase retrieval algorithms: a comparison. Applied optics, 21(15), 2758-2769. (pdf)
• [PB14] Parikh, Neal, and Stephen Boyd. “Proximal algorithms.” Foundations and trends® in Optimization 1.3 (2014): 127-239. (pdf)
• [KZG21] Kirkove, M., Zhao, Y., & Georges, M. P. (2021, July). Inverse-problem-based algorithm for sparse reconstruction of Terahertz off-axis holograms. In Digital Holography and Three-Dimensional Imaging (pp. DM1B-4). Optical Society of America. (url)
• [KZG21b] M. Kirkove, Y. Zhao, M. Georges, “Inverse-problem-based algorithm for sparse reconstruction of Terahertz off-axis holograms"In Digital Holography and Three-Dimensional Imaging (pp. DM1B-4). Optical Society of America. (slides)
• [Baraniuk07] Baraniuk, R. G. (2007). Compressive sensing [lecture notes]. IEEE signal processing magazine, 24(4), 118-121. (pdf)
• [JDCP11] Jacques, L., Duval, L., Chaux, C., & Peyré, G. (2011). A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity. Signal Processing, 91(12), 2699-2730. (pdf)