Image Reconstruction of an Emerging Optical Imager

Abstract

The emerging optical imager can reduce the system weight, size and power consumption by an order of magnitude compared to conventional optical telescopes at the same resolution. It utilizes Fourier-domain interferometry and samples by radial scheme. To date, there is little research on efficient image reconstruction algorithm of it. In this paper, we analyze the sampling properties of the imager and propose a compressive sensing (CS) based image reconstruction method for it. We also discuss the accurate and fast transform from the radial sampling scheme to Cartesian coordinate. The simulations of the proposed method have achieved good performance.

1 Introduction

In the field of astronomy and remote sensing, it is important to obtain the high resolution image of remote target. However, a high resolution telescope requires a large aperture, which comes with a large size, heavy weight and high power. In addition, launching heavy and large objects into space is expensive. These factors restrict the development of large aperture telescope in aerospace. The emerging imaging technology called compact passive coherent imaging technology with high resolution (CPCIT) provides a way to solve this problem.

CPCIT utilizes interference and photonic integrated circuit (PIC) technology to imaging. It can reduce the system weight, size and power consumption by an order of magnitude compared to conventional optical telescopes at the same resolution. In 2010, Fontaine et al. [18] outlined the possibility of amplitude- and phase-accurate measurements. In 2013, researchers at UC Davis and LM Advanced Technology Center developed the segmented planar imaging detector for EO reconnaissance (SPIDER) [13,14,15, 20]. The research was developed by funding from National Aeronautics and Space Administration (NASA) and applied this technology for the first time. To date, only SPIDER has made imaging verification tests [3]. SPIDER is made up of many lenslets arranged in a radial-spoke pattern, each spoke is a 1D lenslets array composing of a set of lenslets. These lenslets comprise different baselines (lenslet pair) and couple light into PIC. PIC is behind the lenslets array and contains optical components for complex visibility (amplitude and phase) measurement. These optical components include arrayed waveguide gratings (AWG), spectral separation, beam combination and $90^\circ$ optical hybrids, etc. Various lenslet pairs are equivalent to many small Michelson stellar interferometers. The SPIDER samples a target that is imaged in the Fourier domain and uses an inverse Fourier transform to reconstruct an image. The schematic diagram of the lenslets arrangement of SPIDER is illustrated in Fig. 1(a). Figure 1(b) and (c) are the first and second generation schematics of PIC. An example set of SPIDER parameters is shown in Table 1.

Fig. 1.

The SPIDER imager concept. Figures are taken from [3]. (a) Lenslet array arrangement of SPIDER. (b) Schematic view of first generation PIC. (c) Schematic view of second generation PIC

As discussed above, CPCIT samples in Fourier domain and requires image reconstruction. To date, there is little research work on image reconstruction for CPCIT. In the paper, we analyze the sampling structure and characteristics of CPCIT and propose an image reconstruction algorithm suitable for CPCIT. The paper is organized as follows. Section 2 discusses the sampling properties of CPCIT and describes the emerging direct and fast Polar Fourier transform, which is suitable for CPCIT. Section 3 presents the compressive sensing based image reconstruction method. The experimental simulations for verifying the proposed method are shown in Sect. 4. Finally we provide a summary and a discussion in Sect. 5.

2 Radial Sampling Scheme

The characteristics of imaging system are greatly influenced by point spread function (PSF), which is closely related to the sampling trajectory. Therefore, studying the sampling scheme of imaging system is very important for image reconstruction. As shown in Fig. 1, CPCIT is sampling in radial scheme, which suits to be described in Polar coordinate system. Besides, there are Cartesian scheme (e.g., most of cameras), spiral scheme (e.g. radio interferometry), etc. Figure 2 shows a Cartesian grid and a Polar grid. The general CPCIT equation for imaging reads as:

$$\begin{aligned} K(\mathbf {\mathbf {u}})=\int _{\mathbb {R}^2}I(\mathbf {\mathbf {a}})exp(i2\pi \mathbf {\mathbf {a}}\cdot \mathbf {\mathbf {u}}) d\mathbf {\mathbf {a}}=\mathcal {F}{\{I(\mathbf {\mathbf {a}})\}}. \end{aligned}$$

(1)

Fig. 2.

Illustration of (a) the Cartesian sampling scheme and (b) the Polar scheme.

Table 1. Example set of SPIDER parameters [16].

Here $\mathcal {F}{\{\cdot \}}$ is the usual Fourier transform. $K(\mathbf {\mathbf {u}})$ is the complex visibility measurement. $I(\mathbf {\mathbf {a}})$ is the intensity distribution of an incoherent optical source, $\mathbf {\mathbf {a}}=(a_1,a_2)\in R^2$ denotes a location lies in the optical source region. $\mathbf {\mathbf {u}}=(u,v)\in R^2$ are the corresponding Fourier coordinates defined by $\mathbf {\mathbf {u}}=\mathbf {\mathbf {b}}/\lambda$, where $\mathbf {\mathbf {b}}$ is the baseline separating the two lenslets. Writing the frequency $\mathbf {\mathbf {u}}=(u,v)$ in Polar coordinates as

$$\begin{aligned} u=\rho cos(\theta ), v=\rho sin(\theta ), \end{aligned}$$

(2)

where $\rho$ represents the radius and $\theta$ is the angle. Then Eq. 1 can be represented as

$$\begin{aligned} K(\rho ,\theta )=\mathcal {PF}{\{I(\mathbf {\mathbf {a}})\}}. \end{aligned}$$

(3)

Here $K(\rho ,\theta )$ is defined with Polar variables and $\mathcal {PF}{\{\cdot \}}$ presents the Polar Fourier transform.

There are some unique properties of radial trajectory that is worth paying attention to [7].

Firstly, when ignoring other effects for simplicity, the PSF can be obtained by setting all sample values to one and zero otherwise and processing the data with the usual image reconstruction procedure. In radial trajectory, the center of the PSF is basically unaffected by changes of the radial count. Most target information remains visible even for a significant amount of streaking artifacts. For this reason, radial sampling offers an attractive under-sampling behavior. In contrast, in Cartesian sampling a reduction of the Fourier domain rows leads to either a resolution decline or the occurrence of aliasing effects that render the image useless.

Secondly, each ray of the radial trajectory carries an equal amount of low and high spatial frequencies. In contrast, for Cartesian sampling the important low spatial frequencies are concentrated in a small number of rows, which creates a difficulty in balancing the acquisition of low and high spatial frequencies in under-sampling scheme. In addition, the aliasing artifacts of a radial under-sampling are scattered to different directions. Thus the incoherence of its PSF masks the artifacts inconspicuous and more noise-like. As all radials pass the center of Fourier domain, the number of sample points is much higher for the low frequencies than for high frequencies. This is essentially different to the Cartesian scheme, where all frequencies are equally covered. Although the radial sample distribution increased complications for the image reconstruction, it is actually a reasonable strategy for typical objects that sampling the low frequencies more densely. In fact, many natural images are characterized by an energy concentration in the center region of the Fourier domain.

We can use these properties to study image reconstruction algorithms that is suitable for CPCIT. The algorithms are detailed in Sect. 3. Here we discuss the first step of image reconstruction, transforming the samples from Polar coordinate to Cartesian coordinate. The traditional method from Polar coordinate to Cartesian coordinate is interpolation. Typically, it is implemented either in the frequency domain (e.g. the regridding method [21]) or in the image domain (e.g. the filtered back projection method [23]). The regridding method interpolates the radial data onto a Cartesian grid by convolution in the Fourier domain, and the approximated Cartesian data can be processed by the standard fast Fourier transform afterwards. However, regridding methods involve significant computation. In order to obtain accurate and fast discrete Fourier transform for Polar grids, we use the direct Polar Fourier transform (PFT) in CPCIT, As discussed in literature [1]. The PFT can accomplish direct and exact computational procedure for the true Polar FT (2D). It is fast, devoid of any oversampling (e.g. PPFFT [2]) or interpolations (e.g. USFFT [17]), and no design on parameters (e.g. MLFFT [24]).

The algorithm of PFT is shown as follows. Consider an image g(r, c) of size $(N+1)\times (N+1)$, where N is even, and r, c are row and column indices respectively. The Polar grid of frequencies ${u_x(p,q),v_y(p,q)}$ is defined in the circle inscribed in the fundamental region, $\{u,v\}\in [N+1,N+1)^2$. The Polar Fourier transform is defined by

$$\begin{aligned} G(p,q)=\sum _{c=-N/2}^{N/2}\sum _{r=-N/2}^{N/2}g(r,c)exp(-i\dfrac{2\pi q}{N+1}(rcos\theta _p+csin\theta _q)), \end{aligned}$$

(4)

where $-N/2\le q\le N/2$ represent the points on a radial ray, $\theta _p=p\varDelta \theta ,0\le p\le P-1$ are the angles. $\varDelta \theta =180^\circ /P$ is the angular spacing and P is the total number of radial rays. PFT uses 1D fractional Fourier Transform (FrFT) for the computation of Fourier coefficients on radial rays of the Polar grid [4]. 1D FrFT, also known as the Chirp-Z transform, is given by

$$\begin{aligned} G^\eta (k)=\sum _{q=-N/2}^{N/2}g(q)exp(-i\dfrac{2\pi k\eta q}{N+1}), -N/2\le k\le N/2, \end{aligned}$$

(5)

where $\eta$ is an arbitrary scalar value, g(q) is a discrete 1D signal, $-N/2\le q\le N/2$ and N is an even integer. Here we require few definitions. The 1D Fourier domain scaling of an image along the X axis can be computed by a 1D FrFT along each row using

$$\begin{aligned} G_x^\eta (r,c)=\sum _{q=-N/2}^{N/2}g(r,q)exp(-i\dfrac{2\pi c\eta q}{N+1}), \end{aligned}$$

(6)

and the Fourier domain scaling along the Y axis is computed along each column using

$$\begin{aligned} G_y^\eta (r,c)=\sum _{q=-N/2}^{N/2}g(q,c)exp(-i\dfrac{2\pi r\eta q}{N+1}). \end{aligned}$$

(7)

Here $-N/2\le r,c \le N/2$.

Consider the level one scaling with $\eta =cos(\varDelta \theta )$, which is composed of two solutions achieved using Eqs. 6 and 7. We take another 1D FrFT which with a variable scale for each column to compute the 2D Fourier domain radial date on the X axis scale grid. For the column c, we scale it by a factor $|c|\xi$ where $\xi =sin(\varDelta \theta )$,

$$\begin{aligned} G_{x}^{\eta ,\xi } (r,c)=\sum _{q=-N/2}^{N/2}g_x^\eta (q,c)exp(-i\dfrac{2\pi r|c| \xi q}{N+1}). \end{aligned}$$

(8)

Similarly, for each row r of the Y axis scaled grid situation, we use the same scale factor $|r|\xi$,

$$\begin{aligned} G_{y}^{\eta ,\xi } (r,c)=\sum _{q=-N/2}^{N/2}g_x^\eta (r,q)exp(-i\dfrac{2\pi c|r|\xi q}{N+1}). \end{aligned}$$

(9)

Here $-N/2\le r,c \le N/2$.

We can observe that only the desired points which accurately overlap with the radial points on the Polar grid rays need to be evaluated. There are four symmetrical radial rays, two at an angle $\varDelta \theta$ and $180^\circ -\varDelta \theta$ from the X axis scaled grid (or basically the horizontal lines), and the other two at angles $90^\circ \pm \varDelta \theta$ from the Y axis scaled grid (or the basically vertical lines). The orthogonal rays at $0^\circ$ and $90^\circ$ are not included. Then the complete level one solution of the Fourier data on four radial rays is given: $G(1,q),G(P-1,q),G(P/2-1,q)$ and $G(P/2+1,q)$, which is denoted by $G_1$. For a Polar grid with P is even, the different radial rays oriented in X axis and Y axis are given by $\theta _l=l\times \varDelta \theta , l\ge 0$. The corresponding scale factors are $\eta _l=cos \theta _l$ and $\xi _l=sin \theta _l$. Consequently, the solutions of all levels can be written as

$$\begin{aligned} G_{x}^{\eta _l,\xi _l} (r,c)=\sum _{q=-N/2}^{N/2}g_x^{\eta _l} (q,c)exp(-i\dfrac{2\pi r|c|\xi q}{l(N+1)}), \end{aligned}$$

(10)

and

$$\begin{aligned} G_{y}^{\eta _l,\xi _l} (r,c)=\sum _{q=-N/2}^{N/2}g_y^{\eta _l} (r,q)exp(-i\dfrac{2\pi c|r|\xi q}{l(N+1)}), \end{aligned}$$

(11)

where $-N/2\le r,c \le N/2$ and $1\le l\le L$. The set of scales $S=[\eta _l, \xi _l; 1\le l \le L]$ has the cardinality of 2L. The complete Polar Fourier transform is given by

$$\begin{aligned} G=\bigcup ^L_{l=1}G_l. \end{aligned}$$

(12)

Here $G_l$ are the solutions for different levels l, $1\le l\le L$.

There is no direct inverse Polar Fourier transform. The proposed PFT can be written as a linear operator $T_{PF}$, then the inverse transform can be approached by solving an optimization problem

$$\begin{aligned} \min _x \parallel T_{PF}x-y\parallel ^2_2, \end{aligned}$$

(13)

where x is real vector of size $(N+1)^2\times 1$ representing the signal, y is a complex valued vector of size $P(N+1)\times 1$ representing the transformed signal. $T_{PF}$ is a complex valued matrix of size $P(N+1)\times (N+1)^2$.

As the local sampling density of points in the Polar grid is varied, weights $\varOmega$ are introduced to compensate for its effects. Then we have

$$\begin{aligned} \min _x \parallel \varOmega \circ T_{PF}x-\varOmega \circ y\parallel ^2_2, \end{aligned}$$

(14)

and we solve it iteratively using the relation of least squares (LS) solution

$$\begin{aligned} x_{j+1}=x_j-\mu T_{PF}^H \varOmega ^2\circ (T_px-y), \end{aligned}$$

(15)

where $\mu$ is the step size and j is the iteration number. $T_{PF}^H$ is the adjoin operator of $T_{PF}$. For details, see [1].

In this section, we introduce a direct discrete Polar Fourier transform, which is suitable for CPCIT. AS CPCIT samples in Polar scheme in spatial frequency domain, it is necessary to interpolate into Cartesian coordinate and then inverse Fourier transform, or directly Polar Fourier transform. Because interpolation may cause inaccuracy and high complexity, the PFT introduced in this section is a good choice.

3 Image Reconstruction

As discussed in Sect. 2, the aliasing artifacts of a radial under-sampling are scattered to different directions. Thus the incoherence of its PSF masks the artifacts inconspicuous and more noise-like. In addition, each ray of the radial trajectory carries an equal amount of low and high spatial frequencies. These characteristics of radial sampling are very suitable for CS reconstruction. In this section, we propose a CS-based image reconstruction method for CPCIT.

3.1 Compressive Sensing

The basic principle of CS is that the compressible or sparse signal can be reconstructed from a small amount of measurements, which can even fewer than the number of the Nyquist theory required [11, 12]. CS has important applications in many areas, such as radio interferometry [28], medical imaging [22], etc. Consider a finite signal $x\in R^N$, a real basis $\varPsi \in R^{N\times T}$, the signal can be decomposed as

$$\begin{aligned} x=\varPsi \alpha , \end{aligned}$$

(16)

where $\alpha$ is a sparse representation of x in $\varPsi$. $\varPsi$ could be wavelet transform or Fourier transform, etc. The linear measurement of x can be expressed as

$$\begin{aligned} y=\varPhi x=\varGamma \alpha , \end{aligned}$$

(17)

where $\varPhi \in R^{M\times N}$ represents measurement matrix, $\varGamma \in R^{M\times N}$ represents sensing matrix, $y\in R^{M}$ represents the observation vector. If $M\ll N$, then the system is highly ill-posed. In the framework of CS theory, x can be exactly reconstructed with a high probability when G satisfies the restricted isometry property (RIP) [5]. Various basises satisfy RIP condition, such as Gaussian random basis [9], Toplit basis [8], Bernoulli basis [10], etc.

Fig. 3.

The test images used in simulations and the sampling arrangement of CPCIT. (a) 1951 USAF. (b) Brid’s Nest. (c) 2D arrangement of lenslet array. (d) corresponding spatial frequency coverage.

From the CS perspective, the ill-posed problem can be solved by minimization the $L_1$-norm of the sparse coefficient $\alpha$

$$\begin{aligned} \min _x\parallel \alpha \parallel _1 \ such \; that \ \parallel y-\varGamma \alpha \parallel _2<\epsilon , \end{aligned}$$

(18)

where the $L_1$-norm is given by $\parallel \alpha \parallel _1=\sum _i\mid \alpha _i\mid$, while the $L_2$-norm is given by $\parallel \alpha \parallel _2=(\sum _i\mid \alpha _i\mid ^2)^{1/2}$. The constraint $\epsilon$ is related to the residual noise level.

3.2 CS-Based Method in CPCIT

The linear representation of the image reconstruction process of CPCIT can be expressed as

$$\begin{aligned} y=Ax+n, \end{aligned}$$

(19)

with

$$\begin{aligned} A=CF, \end{aligned}$$

(20)

where x is the original image with size $N=N^{1/2}\times N^{1/2}$, $y\in R^{M}$ is the M visibility samples. $A\in R^{M\times N}$ represents sensing matrix, $F\in R^{N\times N}$ represents the Fourier transform, which could be fast FT or PFT. $C\in R^{M\times N}$ is the rectangular binary matrix implementing the mask characterizing visibility coverage. The vector n represents noise corrupting the measurements. The Total Variation (TV) minimization problem can also be applied because that many natural signals are also compressible or sparse in the magnitude of their gradient. Replacing the $L_1$-norm sparsity prior with the TV-norm of the signal itself

$$\begin{aligned} \min _x\parallel x \parallel _{TV} \ such \; that \ \parallel y-CFx \parallel _2<\epsilon . \end{aligned}$$

(21)

Here the TV-norm is defined by the $L_1$-norm of the gradient of the signal $\parallel x \parallel _{TV}\,=\,\parallel \nabla x \parallel _1$. We also could pose an enhanced $TV{+}$ problem with the prior knowledge of the positivity of the signal as

$$\begin{aligned} \min _x\parallel x \parallel _{TV} \ such \; that \ \parallel y-CFx \parallel _2<\epsilon \ and \ x\ge 0. \end{aligned}$$

(22)

The basic algorithms in radio telescope research, such as CLEAN [19], MEM [25], can also be used for image reconstruction of CPCIT. But the performances of these algorithms are inferior to CS based algorithms [28], so they are no longer detailed here. Compared to the spiral sampling of radio telescope, the radial sampling of CPCIT is more suitable for CS algorithm. It is worth noting that there are many other recovery algorithms, such as orthogonal matching pursuit (OMP) [26], fast iterative shrinkage-thresholding (FIST) [6], etc. We have only discussed several common ones.

Fig. 4.

The reconstructed images of USAF and Nest using various methods. (a) The reconstructed images (from top to bottom) of USAF using direct FT, BP, TV+, respectively. (b) The reconstructed images (from top to bottom) of Nest using direct FT, BP, TV+, respectively.

4 Experimental Simulation

In this section, we simulate the reconstructed images based on direct Fourier transform and the proposed CS methods. One of the test images is the 1951 USAF resolution test chart (called USAF for simplicity) and the other one is the remote sensing image of the Beijing National Stadium, also known as the Bird’s Nest (called Nest for simplicity). They are shown in Fig. 3(a) and (b), respectively. Both of them have same size of $128\,\times \,128$ pixels. Figure 3(c) is the 2D arrangement of lenslet array and Fig. 3(d) is the corresponding spatial frequency coverage. We use the system parameters in Table 1 for the simulations. The radial number is 37. The peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) are introduced as performance metrics to quantitatively evaluate the reconstructed image quality of these methods [27].

Figure 4 shows the reconstructed images using direct FFT and CS methods which include basis pursuit (BP) (Eq. 18) and TV+ (Eq. 22) approaches. The TV+ method has been shown to achieve good reconstruction quality. Figure 5 shows the reconstruction performance of different approaches with various radial numbers. Figure 6 shows the reconstruction performance of different approaches with various sampling number on a single radial. It can be seen that the reconstruction quality become better when the radial number and the sampling number of a single radial are increasing and the TV+ method achieves a good performance. As the number of sample points is much higher for the low frequencies than for high frequencies in radial sampling scheme, the performance of Direct FFT is acceptable and avoiding the influence of the aliasing effect.

Fig. 5.

The performances of different methods with radial number changed on the (a) USAF and (b) Nest.

Fig. 6.

The performances of different methods with sampling number on a single radial changed on the (a) USAF and (b) Nest.

5 Discussion and Conclusions

In this paper, we present a new imaging technology named CPCIT. We analyze the sampling properties of it and then propose a CS based image reconstruction method. CPCIT can reduce the system weight, size and power consumption by an order of magnitude compared to conventional optical telescopes at the same resolution. As sampling scheme is an important part of image reconstruction, we discuss the radial sampling scheme of CPCIT and use PFT to transform the radial Fourier samples of CPCIT into Cartesian coordinate. However, due to the computation of PFT is rather complicated, we have not succeeded in combining it with CS. The combination of PFT and CS is an important part of future work. In addition, we also consider processing the samples in Polar coordinate, and transforming them to Cartesian coordinates at the end. CPCIT is still in the initial stage of research at present. There is little study of image reconstruction algorithm for it. However, the mature application of CPCIT will have a remarkable impact on the development of astronomical observation and remote sensing.