Sparse Image Reconstruction for the SPIDER Optical Interferometric Telescope

TL;DR

This paper explores the application of sparse radio interferometric image reconstruction algorithms to the SPIDER optical interferometric telescope, aiming to enhance high-fidelity imaging for space science.

Contribution

It adapts and demonstrates radio interferometric sparse reconstruction algorithms for optical imaging with the SPIDER telescope, a novel approach in this context.

Findings

Algorithms successfully reconstruct optical sky images from simulated SPIDER data.

Sparse reconstruction improves image quality and fidelity in optical interferometry.

The approach supports the potential of SPIDER for high-resolution space imaging.

Abstract

The concept of a recently proposed small-scale interferometric optical imaging device, an instrument known as the Segmented Planar Imaging Detector for Electro-optical Reconnaissance (SPIDER), is of great interest for its possible applications in astronomy and space science. Due to low weight, low power consumption, and high resolution, the SPIDER telescope could replace the large space telescopes that exist today. Unlike traditional optical interferometry the SPIDER accurately retrieves both phase and amplitude information, making the measurement process analogous to a radio interferometer. State of the art sparse radio interferometric image reconstruction techniques have been gaining traction in radio astronomy and reconstruct accurate images of the radio sky. In this work we describe algorithms from radio interferometric imaging and sparse image reconstruction and demonstrate their application to the SPIDER concept telescope through simulated observation and reconstruction of the optical sky. Such algorithms are important for providing high fidelity images from SPIDER observations, helping to power the SPIDER concept for scientific and astronomical analysis.

Figures

Tables

Table 1. Table 1: SPIDER configuration parameters adopted from [ Kendrick et al., 2013 ] .

Parameter	Value
Spectral Coverage	500-900 nm
Lenslet Diameter	8.75 mm
Longest Baseline	0.5 m
Number of Lenslets per PIC spoke	24
Number of PIC spoke	37
Number of Spectral Bins	10
FoV at 500 (900) nm	35′′ (65′′)
Maximum Resolution at 500 (900) nm	0.7′′ (1.2′′)
Total Measurements	4440

Equations

$$y(u,v)={\int }_{-\infty }^{\infty }x(l,m)a(l,m){\mathrm{e}}^{-2\pi i(lu+mv)}\mathrm{d}l\mathrm{d}m,$$

$$\mathbf{y}=\mathbf{\Phi }\mathbf{x}+\mathbf{n},$$

$$\mathbf{\Phi }=\mathsf{W}\mathsf{G}\mathsf{F}\mathsf{Z}\mathsf{S},$$

$${\mathsf{G}}_{i,\{k_i+j{\}}_K}=d(u_i-(k_i+j)),$$

$${\mathsf{S}}_{i,i}=s\left(\frac{i}{K}-\frac{1}{2}\right),$$

$$p\left(\bm{x}|\bm{y}\right)\propto\exp\bigl{[}-\|\bm{y}-\mathsf{\bm{\Phi}}\bm{x}\|_{\ell_{2}}^{2}/2\sigma^{2}\bigr{]}\exp\bigl{[}-\gamma\|\bm{\mathsf{\Psi}}^{\dagger}\bm{x}\|_{\ell_{1}}\bigr{]}\,,$$

$$\arg \underset{\mathbf{x}}{\max }p\left(\mathbf{x}|\mathbf{y}\right)$$

$$\arg \underset{\mathbf{x}}{\min }$$

$$\arg \underset{\mathbf{x}}{\min }\Vert {\mathsf{\Psi }}^{†}\mathbf{x}{\Vert }_{{\ell }_1}+{\iota }_{{\mathcal{B}}^{ϵ}(\mathbf{y})}(\mathbf{\Phi }\mathbf{x})+{\iota }_{{\mathbb{R}}_{+}^N}(\mathbf{x}),$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^N}{\min }f(\mathbf{\Phi }\mathbf{x})+g(\mathbf{x}).$$

$${\mathbf{x}}^{(k)}$$

$${\mathbf{v}}^{(k+1)}$$

$${\mathbf{z}}^{(k+1)}$$

$${\sigma }_i=\frac{\Vert \mathbf{\Phi }{\mathbf{x}}_{\mathrm{GroundTruth}}{\Vert }_{{\ell }_2}}{\sqrt{M}}\times 10^{-\frac{\mathrm{ISNR}}{20}}.$$

Full text

1 INTRODUCTION

The Hubble Space Telescope (HST) has changed the way astronomers have looked at the Universe. The number of astronomical studies that have used observations from the HST make it one of the most important observatories in history. More than 15,000 articles that have used HST data, in total collecting 738,000 citations.111See https://www.nasa.gov/mission˙pages/hubble/story/index.html However, telescopes such as the HST and its scientific successor, the James Webb Space Telescope (JWST), are extremely heavy and large, while being expensive in cost and power consumption. Nevertheless such next generation optical telescopes like JWST are critical to address astronomy and cosmology science goals such as answering questions about dark matter through weak lensing and understanding the history and formation of our universe.

Recently, the concept of an instrument known as the Segmented Planar Imaging Detector for Electro-optical Reconnaissance (SPIDER) has been developed [Kendrick et al., 2013, Duncan et al., 2015]. The SPIDER is a small-scale interferometric optical imaging device that first uses a lenslet array to measure multiple interferometer baselines, then uses photonic integrated circuits (PICs) to miniaturize the measurement acquisition. The goal of the SPIDER is to reduce the weight, cost, and power consumption of optical telescopes. Furthermore, additional designs have been proposed that could increase the efficiency of imaging using fewer measurements [Liu et al., 2017, Liu et al., 2018]. Recent visibility measurements using lenslet arrays and PICs have shown to match theoretical predictions [Su et al., 2017]. Unlike traditional optical interferometry, the SPIDER telescope can accurately retrieve both phase and amplitude information [Su et al., 2017], making the measurement process analogous to a radio interferometer. Accurate interferometric image reconstruction methods from radio astronomy can thus be applied to image SPIDER observations.

Radio astronomy has a long history of using interferometry to push beyond the limits of resolution and size, at the computational cost of image reconstruction [Ryle and Hewish, 1960]. An interferometer is a device that measures the cross-correlation function of the signals. Interferometric imaging in the radio has proven to be a popular approach between 50 MHz and 100 GHz, with telescopes such as the Very Large Array (VLA) that have antenna arrays spread over 36 kilometers [Thompson et al., 2008]. The cross-correlation between voltages from each pair of antenna is computed to interfere the complex valued measurements known as visibilities. A visibility represents a Fourier coefficient for the sky brightness, with the Fourier coordinate determined by the antenna pair separation. Typically an antenna pair is known as a baseline, with the baseline length corresponding to the antenna separation [Thompson et al., 2008].

Recently, sparse image reconstruction algorithms that exploit developments from the field of convex optimization have shown to improve the quality of reconstructed observations from radio interferometers considerably, on both simulations and real data [Pratley et al., 2018, Cai et al., 2018b, Pratley et al., 2019a, Pratley et al., 2019b, Pratley, 2019]. In this article we take recent developments from radio interferometric imaging and sparse image reconstruction, and put them into the context of the proposed SPIDER instrument. Such methodology would prove useful in future space based telescopes and space missions based on the SPIDER technology (e.g. aerial observations of planetary surfaces). Ultimately it is evident that recent algorithmic developments for radio interferometric imaging can be directly applied to the SPIDER optical interferometer.

In Section 2 we introduce the background and current developments behind the SPIDER concept. Section 3 discusses the calculation of the SPIDER measurement equation. Section 4 introduces the sparse image reconstruction formalism. Section 5 shows image reconstruction from a simulated SPIDER observation.

2 SPIDER

Key to the concept design of the SPIDER is the use of lenslets to collect signals from incoming light. These signals are combined using a PIC to produce an interferometric measurement (visibility), i.e. a Fourier coefficient of the observation. The Fourier coordinates, $(u,v)$, are determined by the separation size in wavelengths (baseline length) between the lenslets that were used to generate the measurement, with larger separations resulting in higher resolution measurements. However, unlike radio interferometry where all possible pairs of antennas in an array can be combined in an observation, lenslets can only be paired once. If there are $N_{l}$ lenslets, the lenslet array will produce $N_{l}/2$ correlations. This differs to the $N(N-1)/2$ correlations expected from a radio array [Thompson et al., 2008, Liu et al., 2018]. To compensate for this lenslets can be combined with the PIC to split the signal into spectral bins (channels), allowing for increased sampling coverage due to variation of baseline length over wavelength. This strategy has been successful in radio astronomy for decades, and is known as multi-frequency synthesis [Thompson et al., 2008].

The concept design of the SPIDER proposed in [Kendrick et al., 2013] is to put a linear array of lenslets onto a PIC card. The PIC cards are mounted as radial spokes on a disc, producing a radial sampling pattern in the $uv$-plane (however, other sampling patterns are considered in [Kendrick et al., 2013]). The proposed operating wavelengths are between 500 nm and 900 nm. The operating wavelength divided by the size of a lenslet (8.75 mm) determines the field of view to be approximately between 0.5 and 1 arc minutes. The longest baseline along a spoke is 0.5 m, which is sensitive to resolutions between 0.65 and 1.2 arcseconds. Parameters of the SPIDER design adopted from [Kendrick et al., 2013] are listed in Table 1, which leads to the $(u,v)$ sampling coverage shown in Figure 1.

3 SPIDER Measurement Model

The SPIDER is an interferometric instrument operating at optical wavelengths. It follows from [Zernike, 1938] that the measurement equation will have the form of a Fourier Transform

[TABLE]

where $(l,m)$ is the coordinates of the image intensity $x$, $a$ is the sensitivity of the instrument over the field of view, and $(u,v)$ is the separation vector between two lenslets and the Fourier coordinate of the measurement [Thompson et al., 2008]. Due to the limited number of lenslets it is not possible to sample the entire $(u,v)$ domain, leading to an ill-posed inverse problem. However, the measurement equation is linear, and we can represent the measurement equation as a linear operator

[TABLE]

with measurements $\bm{y}_{k}=y(\bm{u}_{k},\bm{v}_{k})$ for a sampling pattern $(\bm{u},\bm{v})$ and intensity values $\bm{x}_{i}=x(\bm{l}_{i},\bm{m}_{i})$ for pixels of an image $(\bm{l},\bm{m})$ and $\bm{n}$ is independent and identically distributed (i.i.d.) Gaussian noise.

The measurement operator requires performing a Fourier transform. However, the sampling pattern of SPIDER does not lie on a regular grid so it is not possible to use the speed of the standard Fast Fourier Transform (FFT). Instead, a non-uniform FFT can be used to evaluate the measurement equation, where an FFT is applied, and the visibilities are interpolated off of the FFT grid. This is known as the non-uniform FFT (NUFFT) but is commonly known as degridding in radio astronomy [Fessler and Sutton, 2003, Thompson et al., 2008].

In a discrete setting, let $\bm{x}\in\mathbb{R}^{N}$ and $\bm{y}\in\mathbb{C}^{M}$. Due to the limited field of view for $\bm{x}$, the values of $\bm{y}$ can be estimated using interpolation off the FFT grid. Sinc interpolation suppresses repeating artifacts from outside the imaged region appearing in the model $\bm{x}$. These artifacts are known as aliasing error, this becomes more difficult to suppress for bright sources outside the field of view. However, Sinc interpolation requires convolution of the FFT grid with a Sinc function for each measurement $\bm{y}_{i}$, which requires a large support. In general, it is possible to replace the Sinc function with an interpolation kernel with minimal support such as a Kaiser-Bessel or Gaussian kernel [Fessler and Sutton, 2003], with appropriate apodization correction in the image domain.

A NUFFT can be represented by the following linear operations

[TABLE]

where $\bm{\mathsf{S}}$ represents scaling required due to the attenuation of the interpolation kernel but also scaling due to the sensitivity of the instrument over the field of view, $\bm{\mathsf{Z}}$ is a zero padding operation of the image to up sample the Fourier domain, $\bm{\mathsf{F}}$ is an FFT, $\bm{\mathsf{G}}$ represents the convolution constructed from the interpolation kernel used to interpolate measurements off of the grid, and $\bm{\mathsf{W}}$ are noise weights applied to the measurements, i.e $\bm{\mathsf{W}}_{kk}=1/\sigma_{k}$ for a variance of $\sigma_{k}^{2}$. This measurement operator is described in detail in [Pratley et al., 2018] where it is applied to radio interferometric imaging. $\mathsf{\bm{\Phi}}$ has its adjoint operator $\mathsf{\bm{\Phi}}^{\dagger}$, which consists of applying adjoints of these operators in reverse. In this work we assume that $\bm{y}$ are the weighted measurements and that the dirty map is defined as $\mathsf{\bm{\Phi}}^{\dagger}\bm{y}$.

3.1 Interpolation kernels

The sparse degridding matrix $\bm{\mathsf{G}}$ can be constructed in 1d from a kernel $d(u)$ by

[TABLE]

where $i$ is the index of the measurement $y_{i}$, $k_{i}$ is the closest integer to visibility coordinate $u_{i}-J/2$ (in units of pixels) and is found by the flooring operation, and $j=1\dots J$ are the possible non-zero entries of the kernel. $\{\cdot\}_{K}$ is the modulo-$K$ function, where $K=\alpha\sqrt{N}$ is the dimension of the Fourier grid in 1-D (for notational sake, the 2d Fourier grid is comprised of $K\times K$ samples). For a separable anti-aliasing kernel, it is straight forward to extend the 1d kernel to 2d. The diagonal operator $\bm{\mathsf{S}}$ is calculated in a similar way by

[TABLE]

where $s(x)$ is the reciprocal of the inverse Fourier transform of $d(u)$. In practice, $\bm{\mathsf{S}}$ can be computed analytically if the inverse Fourier transform of the convolution kernel is known.

The Kaiser-Bessel kernel has shown to perform well as an interpolation kernel while using a minimal support size [Fessler and Sutton, 2003, Pratley et al., 2018]. More examples of interpolation kernels are provided in [Jackson et al., 1991] and [Thompson et al., 2008]. For this work we use the Kaiser-Bessel kernel that was used in [Pratley et al., 2018].

4 Sparse Image Reconstruction

The Bayesian statistical inference framework can be used to address the inverse problem (2), as shown in [Cai et al., 2018a]. From Bayes’ theorem, the posterior can then be expressed as

[TABLE]

where $\bm{\mathsf{\Psi}}^{\dagger}\bm{x}$ represents wavelet coefficients, $\|\cdot\|_{\ell_{1}}$ is the sum of absolute values, and $\|\cdot\|_{\ell_{2}}$ is the Euclidean norm. Maximum a posteriori (MAP) estimation is found by choosing the estimate of $\bm{x}$ that will maximize the posterior, which is equivalent to minimizing the negative log posterior, i.e.

[TABLE]

This minimization problem is known as regularized least squares, with the regularization term in the analysis setting being $\|\bm{\mathsf{\Psi}}^{\dagger}\bm{x}\|_{\ell_{1}}$. In many cases $\|\bm{\mathsf{\Psi}}^{\dagger}\bm{x}\|_{\ell_{1}}$ is chosen to penalize the number of parameters that determine $\bm{x}$ and reduce over fitting; moreover, it can also be used to enforce other properties for $\bm{x}$ like smoothness. Furthermore, it is possible to add indicator functions as a prior that can restrict our solution to be real or positive valued, as is done in the constrained problem below. MAP estimation can be solved efficiently using the Forward Backward Splitting algorithm (e.g. [Cai et al., 2018b]).

An issue of using MAP estimation to perform sparse regularization is choosing a proper regularization parameter $\gamma$ (although there are ways to address this; [Pereyra et al., 2015]). The choice of $\gamma$, however, can be avoided after moving from the unconstrained problem in MAP estimation to the constrained problem

[TABLE]

where $\iota$ is the indicator function that restricts $\mathsf{\bm{\Phi}}\bm{x}$ to the set $\mathcal{B}^{\epsilon}(\bm{y})=\{\bm{q}:\|\bm{y}-\bm{q}\|_{\ell_{2}}\leq\epsilon\}$, $\epsilon$ is the error tolerance, and $\iota_{\mathbb{R}_{+}^{N}}$ restricts the solution to be positive.

4.1 ADMM

In the constrained problem the $\ell_{2}$-norm is replaced with an indicator function that restricts the solution to the $\ell_{2}$-ball of radius $\epsilon$. This indicator function is non-differentiable and it is not possible to apply the Forward Backward method to obtain a solution.

Let $f(\mathsf{\bm{\Phi}}\bm{x})=\iota_{\mathcal{B}^{\epsilon}(\bm{y})}(\mathsf{\bm{\Phi}}\bm{x})$ and $g(\bm{x})=\|\bm{\mathsf{\Psi}}^{\dagger}\bm{x}\|_{\ell_{1}}+\iota_{\mathbb{R}_{+}^{N}}(\bm{x})$. We can then consider the optimization problem

[TABLE]

Problem (9) can be addressed by the alternating direction method of multipliers (ADMM) algorithm [Boyd et al., 2011, Yang and Zhang, 2011, Onose et al., 2016]. Starting from an augmented Lagrangian, variables $\bm{x}$ and $\bm{v}$ are minimized alternatively while updating the dual variable $\bm{z}$ (using the dual accent method; [Boyd et al., 2011]) to ensure that the constraint $\bm{v}=\mathsf{\bm{\Phi}}\bm{x}$ is met in the final solution, i.e.

[TABLE]

The above iterations can be evaluated using a combination of proximal operators and the Dual Forward Backward algorithm (allowing the use of more than one wavelet basis). See [Pratley et al., 2019b] for a detailed description of the ADMM algorithm applied here.

5 Reconstructions

In this section we demonstrate reconstruction of simulated SPIDER observations using the ADMM algorithm, where a solution is found from the constrained problem. We use the software package PURIFY222https://github.com/astro-informatics/purify to perform interferometric image reconstruction, powered by the convex optimization package SOPT333https://github.com/astro-informatics/sopt.

To generate the measurement operator used to simulate the observation we use the Kaiser-Bessel kernel with a support size $J=8$ pixels to reduce aliasing error in the ground truth measurements. For reconstruction, we use a measurement operator with a kernel support size of $J=4$ pixels. The number of pixels in $\bm{x}$ are determined by the ground truth image, $\bm{x}_{\rm GroundTruth}\in\mathbb{R}^{N}_{+}$. We do not include the decrease in sensitivity of the SPIDER instrument away from the center of the field, but this can be included in simulations if it is well characterized. To simulate the observation we follow [Pratley et al., 2018] and add i.i.d.  Gaussian noise to the observational data. We define an input signal to noise ratio (ISNR) to determine the standard deviation of the Gaussian noise, where this standard deviation is defined as

[TABLE]

The Fourier sampling pattern of the observation (i.e. the $uv$-coverage) is determined by the design of the SPIDER instrument and the optical spectral coverage. By combining the entire spectra it is possible to increase the sampling coverage, as explained in Section 2. We use the configuration of Table 1 (shown in Figure 1).

The results presented in Figure 2 show that an observation using the proposed SPIDER design can be effectively reconstructed using PURIFY. Reconstruction was performed using a Dirac basis and Daubechies wavelets 1 to 8. While we have used the design from [Kendrick et al., 2013], where the baselines lie on radial spokes, different baseline configurations may lead to higher quality reconstruction. Depending on the structures in the ground truth sky, different baseline configurations will be more effective at sampling the sky, leading to more effective reconstruction of objects and their details. It was recently shown that the theory of compressive sensing might lead to more efficient designs [Liu et al., 2018].

In summary, we adapt recent developments in radio interferometric imaging, leveraging sparsity and convex optimisation, and show that they are effective for imaging SPIDER observations. Moreover, recent developments in efficient uncertainty quantification for radio interferometric imaging can also be adapted for use with SPIDER [Cai et al., 2018b]. The computational performance of these algorithms can be further increased using GPU multi-threading and distribution across nodes of a computing cluster (as implemented in PURIFY already; [Pratley et al., 2019b]).