Mid-frequency MTF compensation of optical sparse aperture system

Chenghao Zhou; Zhile Wang

doi:10.1364/OE.26.006973

1. Introduction

Optical sparse aperture (OSA) imaging systems are capable of producing high-resolution images while maintaining an overall light collection area that is small with respect to a fully filled aperture achieving the same resolution. The application of sparse-aperture imaging techniques to the field of optical remote sensing has been a topic of research for roughly 50 years [1–3]. Such systems have advantages over resolution-equivalent monolithic apertures due to their lower size, weight, volume [3], and cost [4].

The modulation transfer function (MTF) of sparse aperture imaging systems will, in general, have a main lobe of which means mid-band spatial frequencies contain lower energy than that of a fully-filled aperture with the same spatial frequency cutoff [5]. If the mid-frequency MTF becomes too low, the image information at the corresponding spatial frequencies will be deteriorated or lost. Currently there are two main ways to deal with the mid-frequency defect:

1. Aperture optimization method. In designing an OSA system, a good aperture configuration which will cover all the spatial frequencies of MTF is usually designed using optimization algorithm such as genetic algorithm [6]. Many optimization methods have been proposed in previous researches: Stokes and Duncan improved the mid-frequency MTF by changing the radius of some sub-apertures in classical aperture structure [5]. Mugnier and Rousset proposed a method of aperture optimization according to the minimum image reconstruction error criterion [7]; Tcherniavski and Kahrizi optimized the position and size of the sub-aperture by examining the distribution of spatial frequency spectrum in a particular frequency domain [8]; Salvaggio P S and Schott J R proposed a method of optimizing the positions of sub-apertures of OSA system by genetic algorithm [9].
2. Improving the mid-frequency MTF by image restoration method. This kind of method was firstly proposed by Kevin D. Bell and Richard H. Boucher in 1996 [10]. In their paper, the MTF curves before and after compensation are given, but the specific compensation method is not discussed. In 2002, J. R. Fienup compared the restoration effect of several image restoration methods on OSA images, but his paper did not consider the effect of image restoration on MTF [11]. In 2007, Wang Dayong also carried out experimental research on the image restoration of OSA system, but he also did not discuss the effect of image restoration method on enhancing mid-frequency MTF of the OSA system [12]. In 2011, Zhou Zhiwei and Wang Dayong studied the influence of noise on the restoration of optical synthetic aperture images [13]. Although the methods of image restoration can improve the MTF of the system, when the mid-frequency MTF of the sparse aperture system is missing, the missing part of the mid-frequency information cannot be recovered by restoration method.

Previous studies on MTF of OSA system have focused on analyzing the performance of aperture configuration. However, few studies have been carried out on the mid-frequency MTF enhancement combined with image processing technique, which is an indispensable part of modern optical system imaging. Therefore, the mid-frequency MTF compensation of OSA system through image processing needs further research.

In this paper, the mid-frequency MTF compensation method is studied considering two cases: 1. mid-frequency MTF dropping, and 2. mid-frequency MTF missing in OSA imaging system. Considering the first case, an image restoration method is proposed to compensate the mid-frequency MTF of OSA system with large filling factor. For the second case, when the filling factor is too small, the mid-frequency MTF of OSA system is missing. In this case, the missing mid-frequency information cannot be recovered by image restoration method. Therefore, we adopt two optical systems for imaging the mid and high spatial frequency images respectively. Then the two images are combined to compensate the missing mid-frequency information of the OSA system.

2. Principle of mid-frequency MTF decreasing and missing of OSA system

OSA system is used to separate the single aperture of the system and then rearrange them to form sparse aperture. Its imaging principle is still consistent with imaging principle of optical system. The difference between OSA system and single aperture system is that discrete sub-apertures of OSA system sample the wave-front as discrete wave-front, and then the system converges the discrete light wave-front to form diffraction spots (or interference stripes).

According to Abbe imaging principle [14], the coherent imaging process contains two Fourier transform process. From the object to the back focal plane, the diffracted light wave of the object is decomposed into various frequency angular spectrum components, that is the plane wave component propagated in different directions, and the frequency spectrum of the object is obtained on the back focal plane, thus completing the first Fourier transform. From the back focal plane to the image surface, the frequency spectrum is synthesized into the image, so the second Fourier transform is completed.

Because of the limited size of the aperture, not all the frequency components of object pass through the imaging system. The high frequency components which are outside the cut-off frequency will be cut off by the aperture, leading to blurring of the final image and the loss of the details in the object, as shown in Fig. 1 (a). In the OSA imaging system, because of declining of the filling factor, some of the mid-frequency components of objects cannot pass through the imaging system as shown in the Fig. 1(b). The intensity of the image is then weakened after imaging and thus the coherent transfer function (CTF) is decreasing or even missing at the mid-frequency. According to the Fourier optics, for the same optical system, the relationship of optical transfer function (OTF) and CTF can be written as:

O T F (ξ, η) = \frac{C T F (ξ, η) \otimes C T F (ξ, η)}{\iint {| C T F (α, β) |}^{2} d α d β} .

Equation (1) shows that OTF is equal to the normalized autocorrelation of CTF, which is normalized autocorrelation of pupil function [14]. Therefore, discontinuity of the aperture leads to decreasing and missing of mid-frequency OTF in the OSA system.

Fig. 1 Abbe imaging principle: (a) Abbe imaging principle of single aperture system; (b) Abbe imaging principle of OSA system.

Download Full Size | PDF

According to the analysis above, whether the mid-frequency MTF is decreasing or missing is determined by the filling factor of OSA system. When the filling factor is large, the mid-frequency MTF shows a decline, and the high frequency MTF is poor. When the filling factor is small, the mid-frequency MTF is missing, and the high frequency MTF is better. Figure. 2(b) illustrates the relationship between filling factor and MTF in the baseline direction of two sub-apertures system shown in Fig. 2(a). In Fig. 2(b), F denotes the filling factor. The relationship between MTF and filling factor of two sub-apertures system can be expressed as:

M T F = M T F_{d} (f_{x}, f_{y}) + \frac{1}{2} [M T F_{d} (f_{x} + (1 - \sqrt{\frac{F}{2})} ρ_{D c}, f_{y}) + M T F_{d} (f_{x} - (1 - \sqrt{\frac{F}{2})} ρ_{D c}, f_{y})],

where ρ_Dc denotes the cutoff frequency of equivalent aperture. MTF_d is the MTF of sub-aperture which can be described as:

M T F (f) = \frac{2}{π} [arc \cos \frac{f}{f_{c}} - \frac{f}{f_{c}} \sqrt{1 - {(\frac{f}{f_{c}})}^{2}}],

where f_c denotes the cut-off frequency of sub-aperture. Whether there are zero values within the cut-off frequency can be calculated by the feature triangulation method [16].

Fig. 2 Relationship between filling factor and MTF aligned the baseline of the two-aperture system: (a) two sub-aperture system; (b) filling factor versus MTF.

Download Full Size | PDF

3. Mid-frequency MTF compensation for large filling factor OSA system

Now there are two main ways to deal with the mid-frequency defect: one is aperture optimization; another one is image restoration method introduced by Kevin D. Bell and Richard H. Boucherin in 1996. But they did not discuss their method in detail. Figure 3 shows the result of their work. MTFs curves of Golay-6 configuration before and after compensation, and comparison to a monolithic aperture are given in Fig. 3 [10].

Fig. 3 MTF curves before and after compensation [10]

Download Full Size | PDF

The filling factor of the system affects the degradation of mid-frequency MTF of OSA system. When the filling factor is too small, some of mid-frequency MTF of OSA system will drop to zero and the image quality is seriously affected. In this situation, the missing mid-frequency information cannot be recovered by image restoration. So other methods are needed to compensate for the missing mid-frequency information.

According to the analysis in the Ref [10], we can see that post-processing is necessary to restore sparse aperture images, but this process can also give rise to unpleasant artifacts. The main reason for these artifacts is that the MTF (or PSF) used in the image restoration does not match the actual MTF (or PSF). Specifically, the difference between actual and used MTF (or PSF) is mainly caused by wave front error of each FOV.

Because of existence of aberration and difference of aberration in each FOV, the PSF of OSA system is not spatially constant in FOV of the system. The OSA system is a spatial variant system [15]. In addition, errors in the wave front in different FOVs can result in artifacts that degrade the utility of the imagery for either analyzing or testing algorithms [15]. Figure. 4 shows the artifacts caused by the wave front error.

Fig. 4 Restored images: (a) restored image without wave front-error correction; (b) restored image with wave-front error correction [15]. The author has been authorized to use of this figure.

Download Full Size | PDF

Figure 4(a) shows artifacts that arise when wave front error is present but not take into consideration in restoration filter. In addition to the obvious smearing artifacts, the contrast of the image decreased, deteriorating the image quality [15]. Based on the analysis, a method of image restoration based on spatial-variant (SV) PSF is necessary to restore the image of OSA system.

In order to deal with the artifacts caused by the different wave front errors in different FOVs in OSA system, an image restoration method is proposed based on spatial-variant PSF. There are four major problems in spatial-variant PSF based image restoration: (1) the problem of acquiring and storage of SVPSF of each FOV; (2) the problem of fast calculation of deconvolution [21]; (3) the problem of suppressing the noise and ill-posedness of inverse problem [21]; (4) the problem of solving the restoration algebraic equation. These four problems and solutions to them are shown in Fig. 5.

Fig. 5 Framework of SV image restoration.

Download Full Size | PDF

The purpose of this section is to propose a framework (shown in Fig. 5) of a method for realizing SV image restoration algorithms based on SVPSF which will produce some promising results in restoring OSA system.

3.1. Decomposition and acquisition of Space Variant PSF

In direct SV image restoration, all the PSFs of each FOV corresponding to the pixel on image plane are needed. As we know, it is almost impossible to get all the PSFs of each FOV of the optical system. The unknown PSFs can be determined from known PSFs by the means of PSF interpolation and decomposition.

3.1.1. Measurement of PSFs

Because it is hard to obtain all the PSFs of each FOV, we need to measure some PSFs on image plane (shown in Fig. 6(a)) in order to obtain all the PSFs by means of interpolation and decomposition. In practical application when an actual system is built, the known PSFs should be obtained by the means of point pulse method using the target plane shown in the Fig. 6(b). Some PSFs at different FOV can be acquired by point pulse method on image plane as shown in Fig. 6(a). This target plane shows the position of the point pulse on the object plane of OSA imaging system. The target plane should cover the entire FOV of the OSA system. In simulation experiment in this section all the known PSFs are simulated in Zemax.

Fig. 6 Position map of PSF (a): position of measured PSF on image plane; (b): target plane put on object plane.

Download Full Size | PDF

3.1.2. Unknown PSFs acquisition and storage by decomposition and interpolation

After measuring some of the PSFs, we need to obtain all the PSFs on the image plane. The basic idea of obtaining PSFs by interpolation and decomposition is that a set of basis function should be found, which decompose the known PSFs on the image plane and express the known PSFs by these set of basis functions and the corresponding coefficients. Then by interpolation, we can calculate the coefficients of unknown FOV. Then bring the coefficients into the basis functions to acquire the unknown PSFs. So there are two steps in this subsection: PSFs decomposition and coefficients interpolation.

1) PSFs decomposition based on PCA method
A lot of methods can be used to decompose PSFs such as Karhunen-Loeve (K-L) decomposition, wavelets basis functions, Shapelet basis functions, principal components analysis (PCA) and so on. In this paper PCA method is used for decomposing PSFs. Because PCA method is simple and convenient. And it can build the basis functions of PSF according to optical system. So the basis function of PSF of a certain optical system obtained by PCA method can have better adaptation to this optical system. In the step of interpolation, PSFs obtained by PCA method also have very good precision which meets the requirement of image restoration. The PCA method was first proposed by Hotelling in 1933 [17]. The idea of this method is to convert multi-dimensional data into low-dimensional data which retains the main information provided by the original data and not related to each other. By using the PCA method, it is easy to grasp the main features of complex data set. For these reasons, PCA method appears to be an attractive way to decompose the known PSFs. The remainder of this subsection will address the PCA processes in great details.
The process of PSF decomposition using PCA method can be represented as follows.
- (1) Normalize measured PSFs: ${PSF}_{i} = \frac{{PSF}_{i}}{sum ({PSF}_{i})}, i = 1, 2, 3, ... N,$
  where sum denotes the sum of all the elements of PSF_i.
- (2) Rewrite all the known PSFs as one-dimensional vectors and then standardize them by: ${PSF}_{i}^{*} {=PSF}_{i} - mean (PSF), i = 1, 2, 3, ..., N,$
  where mean denotes the mean of all the known PSFs.
- (3) Calculate covariance matrix of known PSFs: $C_{i j} {=cov(PSF}_{i}^{*}, {PSF}_{j}^{*}), 1 \leq i, j \leq N,$
  where cov denotes covariance of PSF_i* and PSF_j*.
- (4) Then eigenvalues {λ₁, λ₂, λ₃, …}, λ_N} and corresponding eigenvectors {eig₁, eig₂, eig₃, …, eig_N} of covariance matrix are calculated.
- (5) In order to select number of basis functions, we need to calculate principal component scores: $ρ_{i} = \frac{λ_{i}}{\sum_{j = 1}^{n} λ_{j}},i=1,2,3, ...,N,$
  where ρ_i denotes cumulative contribution rate. In general, the principal components, which make the cumulative contribution rate more than 85% [17], should be chosen in order to return the original data set and maintain original information.
- (6) By using the eigenvectors of covariance matrix the basis functions {r₁, r₂, r₃, … r_N} are constructed: $(\begin{matrix} r_{1} \\ r_{2} \\ r_{3} \\ ... \\ r_{N} \end{matrix}) = (\begin{matrix} e i g_{1} \\ e i g_{2} \\ e i g_{3} \\ ... \\ e i g_{N} \end{matrix}) (\begin{matrix} {PSF}_{1}^{*} \\ {PSF}_{2}^{*} \\ {PSF}_{3}^{*} \\ ... \\ {PSF}_{N}^{*} \end{matrix}) .$
- (7) PSF data regression using the basis function and coefficients: $(\begin{matrix} {PSF}_{1}^{*} \\ {PSF}_{2}^{*} \\ {PSF}_{3}^{*} \\ ... \\ {PSF}_{N}^{*} \end{matrix}) = {(\begin{matrix} e i g_{1} \\ e i g_{2} \\ e i g_{3} \\ ... \\ e i g_{N} \end{matrix})}^{T} (\begin{matrix} r_{1} \\ r_{2} \\ r_{3} \\ ... \\ r_{N} \end{matrix}) + mean (PSF),$
  where {eig₁, eig₂, eig₃, …, eig_N} is the coefficients matrix of the basis function.
2) Coefficients of unknown PSFs obtained by interpolation
The basic idea of PSF acquiring by interpolation is that PSF at each interpolation point has a certain basis function representation; then according to the coefficients of the basis functions and coordinates of known points, the coefficients of basis functions of unknown points are obtained by interpolation. Thus PSFs of unknown points are constructed by acquired coefficients and basis functions:
$PSF = (\begin{matrix} a_{1} & a_{2} & a_{3} & ... & a_{k} \end{matrix}) (\begin{matrix} r_{1} \\ r_{2} \\ r_{3} \\ ... \\ r_{k} \end{matrix}) + mean (PSF),$
where {a₁, a₂, a₃, …, a_N} is the coefficients obtained by interpolation

Many kinds of interpolation methods can be used to get coefficients of unknown PSFs: inverse distance to a power, Kriging, minimum curvature, modified Shepard’s method, natural neighbor, nearest neighbor, polynomial regression, radial basis function, triangulation with linear interpolation, moving average and local polynomial [18].

By PCA decomposition and interpolation, the PSF of each FOV can be decomposed into a set of basis functions and the corresponding coefficients. So PSFs data of all the FOVs can be stored as a set of basis functions and the corresponding coefficient matrixes, that is to say only several matrixes (often no more than twelve) need to be stored in the computer. This greatly reduces the storage space. We can get PSFs of all FOVs easily.

3.2. The SV image restoration algorithm

3.2.1. Calculation of deconvolution

In order to decrease computational complexity, we need to deal with the calculation of deconvolution. There are several methods to deal with calculation of deconvolution such as circulant matrix model and aperiodic model. In 1971, B. R. Hunt and R. M. Gray published a number of papers which suggested and promoted the circulant matrix model [19–21]. In this method, the imaging process can be written as:

y = B_{c} x,

where B_c denotes the circulant matrix. According to the properties of circulant matrix, B_c can be approximated as:

B_{c} {=WΛW}^{-1},

where W is formed by eigenvector of the circulant matrix; Λ is the Fourier transform of the PSF matrix. And W also has the following properties:

WΛ \Leftrightarrow M N • IDFT[Λ],

W^{-1} Λ \Leftrightarrow \frac{1}{M N} • DFT[Λ],

where M, N denotes the size of the Λ; DFT and IDFT denotes the discrete Fourier transform operation and inverse discrete Fourier transform operation. We can see that convolution process can be understood as the process of inverse Fourier transform and Fourier transform. Image restoration algorithm can be calculated by fast Fourier transform algorithm.

As shown in Eq. (10), SVPSFs can be written as:

PSF(x,y) = \sum_{k = 1}^{{PCA}_{num}} a_{k} (x,y) • r_{k} (x,y) + mean (PSF),

where PCA_num denotes the number of the basis functions. So the image process of SV system can be written as:

\begin{matrix} y = B_{c} x = \sum_{i = 1}^{M N} x_{i} * P S F_{i} \\ = \sum_{i = 1}^{M N} x_{i} * (\sum_{k = 1}^{{PCA}_{num}} a_{k} (x,y) • r_{k} (x,y) + mean (PSF)) \\ = \sum_{i = 1}^{M N} x_{i} * \sum_{k = 1}^{{PCA}_{num}} a_{k} (x,y) • r_{k} (x,y)+ \sum_{i = 1}^{M N} x_{i} * mean (PSF) . \\ = \sum_{k = 1}^{{PCA}_{num}} r_{k} * (A_{k} x) + x *mean(PSF) \end{matrix}

where A_k is the matrix of coefficients a_k and i denotes the ith FOV on the image plane. According to the circulant matrix model:

\begin{matrix} \sum_{k = 1}^{{PCA}_{num}} [A_{k} x * r_{k}]+ x *mean(PSF) \\ = \sum_{k = 1}^{{PCA}_{num}} [A_{k} x *W Λ_{k} W^{- 1}]+ x *W Λ_{m e a n} W^{- 1} . \end{matrix}

So the circulant matrix B_c of the SV image process is:

B_{c} = \sum_{k = 1}^{{PCA}_{num}} [A_{k} W Λ_{k} W^{- 1}]+W Λ_{m e a n} W^{- 1} .

3.2.2. Regularization of deconvolution problem

Due to the influence of noise, image deconvolution is an ill-posed inverse problem. One way to solve the ill-posed inverse problems is by using a process called regularization. This process introduces additional information in the form of extra constraints. Methods such as constrained least squares method, Tikhonov regularization and total variation method can also be good choices to solve ill-posed inverse problem. In this paper, we choose constrained least squares method as an example method of solving ill-posed inverse problem. We just give the restoration equation of the method, for more detailed information, see Ref [21]. The restoration equation is:

(B_{c}^{T} B_{c} + α C^{T} C) x = B_{c}^{T} y,

where C is the Laplace operator:

c (m, n) = \frac{1}{8} [\begin{matrix} 0 & 1 & 0 \\ 1 & - 4 & 1 \\ 0 & 1 & 0 \end{matrix}] .

And α determines the effect of the regularization term. The value of alpha usually between 0 and 1. There are some methods developed to choose the value of alpha such as Generalized Cross Validation, the Discrepancy Principle and the L-curve method [22]. All this method can be used to choose value of alpha in this paper. In the simulation experiment, we can also manually selected the values of alpha based on the image quality evaluated by ISNR to assist α selection:

I S N R = 10 \log_{10} \frac{{‖ g - f ‖}^{2}}{{‖ f - f_{α} ‖}^{2}}

where g denotes the original image, f denotes the blurred image and f_α is the restored image using a certain value of alpha. The most suitable alpha would be the one making the ISNR maximum.

The matrix B_c of Eq. (18) can be taken into the Eq. (19):

\begin{matrix} (\sum_{k = 1}^{{PCA}_{num}} [A_{k} W Λ_{k} W^{- 1}]+W Λ_{m e a n} W^{- 1})^{T} • \\ (\sum_{k = 1}^{{PCA}_{num}} [A_{k} W Λ_{k} W^{- 1}]+W Λ_{m e a n} W^{- 1}) x + α C^{T} C x, \\ = (\sum_{k = 1}^{{PCA}_{num}} [A_{k} W Λ_{k} W^{- 1}]+W Λ_{m e a n} W^{- 1})^{T} y \end{matrix}

Equation (22) is the constrained least squares method of SV image restoration algebraic equation.

3.2.3. Solving of algebraic equation

The solving of algebraic equation is also an important question. The use of regularization technique will lead to a minimization problem, which is usually solved by iteration method. In this paper, we choose conjugate gradient iteration to solve Eq. (22).

3.3. Simulation experiment

In order to compare the recovery ability of space invariant restoration algorithm and space variant restoration algorithm, in this paper we carried out the simulation experiment. In simulation experiments, all the PSFs data and simulation images were obtained from Zemax and Matlab. We first designed a phased array OSA system using Zemax as shown in Fig. 7. It was made of three groups of system: afocal telescope, beam combiner telescope and periscopic active mirror.

Fig. 7 Phased array OSA system: (a) shaded model of OSA system; (b) 3D schematic layout.

Download Full Size | PDF

Three images were selected in this simulation experiment. And according to the PSF acquired in Zemax we simulated nine images whose size are 256 × 256 and signal to noise ratio (SNR) are 20dB, 30dB and 50dB, respectively. The Gaussian white noise are used in this simulation experiment. Real systems will have more structured noise that may have a higher impact on image quality. SNR are calculated by the following equation:

S N R = 10 \log_{10} (\frac{var (g)}{σ_{n}^{2}}),

where var(g) is the variance of blurred image without noise;

σ_{n}^{2}

is the variance of Gaussian white noise. Three simulated images are shown in Fig. 8.

Fig. 8 Simulated space variant images: (a) SNR = 20; (b) SNR = 30; (c) SNR = 50.

Download Full Size | PDF

The framework of our space variant constrained least squares algorithm is shown in Fig. 5. Four kinds of interpolation methods are used to get the coefficients of the unknown PSF values: polynomial regression, inverse distance to a power, modified Shepard’s method, radial basis function [18]. Figure 9 shows the accuracy of the acquired PSF by means of interpolation. The accuracy of each PSF of FOV can be calculated by mean square error (MSE) which can be express as:

M S E (i, j) = {\sum_{k = 1}^{m} \sum_{l = 1}^{n} (P S F_{zemax} (k, l) - P S F_{i p} (k, l))}^{2},

where m, n denotes the size of the PSF; i, j denotes the position of PSF on the image plane, and average MSEave is calculated by:

{MSE}_{a v e} = (\sum_{i = 1}^{M} \sum_{j = 1}^{N} M S E (i, j)) / (M N)

where M, N denotes the size of the image, the average MSE_ave of each interpolation method are shown in Table 1. We can see from Fig. 9. that the MSE of PSFs of edge FOV is lower than that of center FOV in all the interpolation methods. We can conclude that polynomial regression method, modified Shepard’s method and radial basis function method all have a good accuracy The MSE of PSFs acquired by these three methods can reach the order of magnitude of 10⁻⁶. Obviously in Table. 1, the MSE_ave of inverse distance to a power method is worse than other methods. The accuracy of radial basis function method is the best one, the MSE_ave of which is 1.396 × 10⁻⁶. In these simulation experiments the size of the PSF matrix is 31 × 31.

Fig. 9 Accuracy map of PSF acquired by interpolation (a): polynomial regression; (b): inverse distance to a power; (c): modified Shepard’s method; (d): radial basis function.

Download Full Size | PDF

Table 1. Average MSE of each interpolation method.

View Table | View all tables in this article

Table 2 shows the improved signal to noise ratio (ISNR) of each group of restored images. ISNR are calculated as:

Table 2. ISNR results of each group of restored image.

View Table | View all tables in this article

I S N R = 10 \log_{10} \frac{{‖ g - f ‖}^{2}}{{‖ f' - f ‖}^{2}}

From the Table 2, we can conclude that the restored image has better quality and higher ISNR is achieved. In addition, as the ISNR cannot be a perfect index to judge the image quality, we also show the local enlarged image of the restored images. Due to the article space limitations, we only show the part of the simulated image of 30dB and part of restored image of 30dB, as shown in the Fig. 10. We can see both from the Table 2 and from Fig. 10 that the quality of the space variant restored image is better than the space invariant restored image. We can see that the sharpness of the surveillance camera (as shown in the red box) in Fig. 10(d) is better than that in Fig. 10(c).

Fig. 10 Restored SIV images: (a) part of original image; (b) part of blurred image; (c) part of SIV restored image; (d) Part of SV restored image.

Download Full Size | PDF

4. Mid-frequency MTF compensation for small filling factor OSA system

The image restoration method is used to compensate the mid-frequency MTF of OSA system with large filling factor. When the filling factor is too small, the mid-frequency MTF of OSA system is dropping to zero and therefore mid-frequency information are not obtained. In this situation, the missing mid-frequency information cannot be recovered by image restoration. In this paper, two optical systems are used for imaging the mid and high spatial frequency images respectively, and then the two images are fused to compensate the missing mid-frequency information of the OSA system.

4.1 Method of compensating missing mid-frequency information of OSA system

According to the analysis in Sec. 2, when the filling factor is small, the MTF of OSA system exhibits stopband characteristics as shown in Fig. 2, and it can be seen that the smaller the filling factor is, the better the MTF of high-frequency is.

Imaging is the process of obtaining information of different object spatial frequencies, which can be simultaneous or not. Thus, two optical systems are used in this paper for imaging the mid and high spatial frequency images respectively, and then the two images are fused by pixel level image fusion method to compensate the missing mid-frequency information in OSA system image. This method is shown in Fig. 11.

Fig. 11 Compensating missing mid-frequency information of OSA system.

Download Full Size | PDF

4.2 Feasibility of compensation of missing mid-frequency MTF

According to Fourier optics and geometrical optics, the intensity distribution function of object is a periodic function. It can be seen as superposition of a set of cosine function of different frequencies, amplitudes and initial phases. For non-periodic functions, it can be regarded as a periodic function whose period tends to be infinite. That is to say, for one dimensional (two-dimensional images are similar to one dimensional images) discrete image, the imaging process can be expressed as:

g (x) = \sum_{0}^{f_{c u}_{t}} M T F (f) A (f) \cos (f x + θ_{f}),

where g(x) denotes the images, f denotes the spatial frequency, A(f) denotes the amplitudes and θ_f denotes the initial phase. The imaging process of OSA system can be expressed as:

g_{O S A} (x) = \sum_{0}^{f_{low}} M T F_{l o w} (f) A (f) \cos (f x + θ_{f})+ \sum_{0}^{f_{h i g h}} M T F_{h i g h} (f) A (f) \cos (f x + θ_{f}),

where MTF_low and MTF_high are the low-frequency and high-frequency parts of the MTF of the OSA system. The imaging process of complementary single aperture system, which is used to compensate the mid-frequency MTF of OSA system can be expressed as:

g_{s i n g l e} (x) = \sum_{0}^{f} M T F_{s i n g l e} (f) A (f) \cos (f x + θ_{f}),

where MTF_single denotes the MTF of the single aperture system. According to Eq. (28) and Eq. (26), it can be seen that the whole spatial frequency image can be obtained simply by adding the two images or other image fusion methods.

4.3 MTF analysis of OSA system and single aperture system

In this paper, we analyze the MTF of an OSA system with four sub-apertures with filling factor of 0.333 and a single aperture system as an example. The aperture structure of the OSA system is shown in Fig. 12.

Fig. 12 Aperture structure of the OSA system.

Download Full Size | PDF

The filling factor and cut-off frequency of the OSA system are as follows [6]:

F = \frac{4 π d^{2}}{π D^{2}},

f_{D c} = \frac{D}{λ f},

where D and d denotes the diameter of encircled aperture and sub-aperture, λ denotes the wavelength and f denotes the focal length. MTF of the OSA system can be expressed as [6]:

\begin{matrix} M T F (f_{x}, f_{y}) = \frac{1}{4} \sum_{i = 1}^{4} \sum_{j = 1}^{4} M T F_{d} (f_{x} - \frac{x_{i} - x_{j}}{λ f}, f_{y} - \frac{y_{i} - y_{j}}{λ f}) \\ = M T F_{d} (f_{x}, f_{y}) + \frac{1}{4} [M T F_{d} (f_{x} \pm \frac{D - d}{λ f}, f_{y}), \\ + 2 M T F_{d} (f_{x} \pm \frac{D - d}{2 λ f}, f_{y} \pm \frac{D - d}{2 λ f}) \\ + M T F_{d} (f_{x}, f_{y} \pm \frac{D - d}{λ f})] \end{matrix}

where MTF_d denotes the MTF of sub-aperture,

d = \frac{\sqrt{F}}{2} D

and

\frac{(D - d)}{λ f} = (1 - \frac{\sqrt{F}}{2}) f_{D c}

. So Eq. (29) can be simplified as [6]:

\begin{matrix} M T F (f_{x}, f_{y}) = \frac{1}{4} \sum_{i = 1}^{4} \sum_{j = 1}^{4} M T F_{d} (f_{x} - \frac{x_{i} - x_{j}}{λ f}, f_{y} - \frac{y_{i} - y_{j}}{λ f}) \\ = M T F_{d} (f_{x}, f_{y}) + \frac{1}{4} [M T F_{d} (f_{x} \pm (1 - \frac{\sqrt{F}}{2}) f_{D c}, f_{y}) \\ + 2 M T F_{d} (f_{x} \pm \frac{(2 - \sqrt{F}) f_{D c}}{4}, f_{y} \pm \frac{(2 - \sqrt{F}) f_{D c}}{4}), \\ + M T F_{d} (f_{x}, f_{y} \pm (1 - \frac{\sqrt{F}}{2}) f_{D c})] \end{matrix}

According to the Fourier optics, the missing frequency range [f_x₁, f_x₂] of the OSA system in the horizontal and the vertical direction is:

f_{x 1} = \frac{d}{λ f} = \frac{\sqrt{F} D}{2 λ f} = \frac{\sqrt{F}}{2} f_{D c},

f_{x 2} = \frac{(D - 2 d)}{λ f} = (1 - \sqrt{F}) \frac{D}{λ f} = (1 - \sqrt{F}) f_{D c} .

The MTF of single aperture system is as follows:

M T F (f_{x}, f_{y}) = \frac{2}{π} [arc \cos (\frac{\sqrt{f_{x}^{2} + f_{y}^{2}}}{f_{0}}) - (\frac{\sqrt{f_{x}^{2} + f_{y}^{2}}}{f_{0}}) \sqrt{1 - {(\frac{\sqrt{f_{x}^{2} + f_{y}^{2}}}{f_{0}})}^{2}}],

where f₀ = D_x/λf is the spatial cut-off frequency of the optical system, D_x is the diameter of the single aperture system. It is necessary to compensate the missing mid-frequency MTF of the OSA system. The MTF of the missing mid-frequency lower limit f_x₂ should be at least 0.3 (based on optical design experience) or even higher:

M T F (f_{x}, 0) = \frac{2}{π} [arc \cos (\frac{f_{x 2}}{f_{0}}) - (\frac{f_{x 2}}{f_{0}}) \sqrt{1 - {(\frac{f_{x 2}}{f_{0}})}^{2}}] = 0.3.

The solution D_x of the Eq. (33) is the diameter of the single aperture system.

4.4 Image fusion method

According to the analysis in Sec. 4.1, the whole spatial frequency image can be obtained by simply adding the images from the two systems directly. So the weighted image fusion method is used:

F (x, y) = ω_{1} A (x, y) + ω_{2} B (x, y),

where ω₁ and ω₂ is the weighted coefficients, ω₁ + ω₂ = 1.

The purpose of the image fusion is to compensate the mid-frequency MTF of the OSA system. Some real images may not contain the target spatial frequency information, thus the average weighted image fusion method cannot obtain the image of maximum spatial frequency information. Therefore, the weights need to be determined according to the amplitude of the spatial frequency contained in the two images:

a = \int_{f_{x 2}}^{f_{_{D c}}} A (f) d f,

b = \int_{f_{x 1}}^{f_{x 2}} B (f) d f,

ω_{1} = \frac{a}{a + b},

ω_{2} = \frac{b}{a + b},

where A and B denotes the spatial spectral intensity of the two images; a represents the intensity integral of the high frequency spectrum of OSA image, and b represents the integral of the mid-frequency spectrum of single aperture image.

4.5 Simulation experiments

In this section, Zemax and Matlab software are used to test the proposed method. The radial target image is used in this simulation experiment. The radial target clearly illustrates the missing mid-frequency OSA system. The OSA system with filling factor of 0.3333 and the diameter of the equivalent aperture of 540mm was designed in Zemax (shown in Fig. 13). According to Eq. (12), the diameter of the single aperture system is 390mm. In this example, the equivalent aperture is increased about 150mm that is about increasing by 38.5% compared with the single aperture system. It has a large increase and has a certain value of application.

Fig. 13 Optical structure of the two systems in Zemax: (a) OSA system; (b) single aperture system.

Download Full Size | PDF

The PSFs and MTFs (shown in Fig. 14) of two designed systems, are extract in Zemax and used for image simulation. Figure 14 shows the MTF of the two system in the horizontal and vertical direction, as well as the MTF of the equivalent aperture system.

Fig. 14 MTF curves of the two systems and MTF curve of the equivalent aperture system.

Download Full Size | PDF

Figure 15(a) and Fig. 15(b) are the images of OSA system and single aperture system. It can be seen from the figure: in Fig. 15(a), mid-frequency stripes in horizontal and vertical direction are missing (as shown in the red box) while the high stripes are clear; and in Fig. 15(b) the high frequency stripes are lost (as shown in the red box) while mid-frequency stripes exists. Figure 15(c) is the equivalent aperture system image with mid and high frequency stripes in it. Figure 15(d) is the fused image, it can be seen that the fused image contains both the mid frequency and the high frequency stripes in the horizontal and vertical direction.

Fig. 15 Results of simulation experiment: (a) Image of OSA system; (b) Image of single aperture system; (c) Image of equivalent aperture system; (d) fused image.

Download Full Size | PDF

5. Summary and conclusion

Sparse aperture imaging is a useful technique for obtaining high resolution images while maintaining a small light collection area with respect to a large, resolution-equivalent single aperture. However, conventional OSA system suffer from reduced mid-frequency contrast.

In this paper, the principle of mid-frequency MTF decreasing and missing are analyzed. Then, according to the principle of the decreasing and missing of the mid-frequency MTF, two different mid-frequency MTF compensation methods are proposed.

Simulation experiments were performed to examine the methods of compensation missing and decreasing mid-frequency MTF. These numerical experiments verified with little disagreement between theory and simulation results. As the results will compensate the missing and decreasing mid-frequency MTF of OSA system, we believe that our mid-frequency MTF compensation method will be useful to other researchers examining a wide variety of other OSA system.

Funding

National Natural Science Foundation of China (NSFC) (61007008);

Acknowledgments

We gratefully acknowledge Changliang Guo, Jiajie Wu and Min Lu. And we also acknowledge further support from the Research Center for Space Optical Engineering of Harbin institution of technology.

References and links

1. A. B. Meinel, “Aperture Synthesis Using Independent Telescopes,” Appl. Opt. 9(11), 2501 (1970). [CrossRef] [PubMed]

2. S.-J. Chung, D. W. Miller, and O. L. de Weck, “Design and implementation of sparse aperture imaging systems,” Proceedings of SPIE - The International Society for Optical Engineering 4849, 181–192 2002 (2002). [CrossRef]

3. J. S. Fender, “Synthetic apertures: An Overview,” Proceedings of SPIE - The International Society for Optical Engineering 440, 2–7 (1984). [CrossRef]

4. J. R. Fienup, “MTF and integration time versus fill factor for sparse-aperture imaging systems,” Proc. SPIE 4091, 43–47 (2000). [CrossRef]

5. A. J. Stokes, B. D. Duncan, and M. P. Dierking, “Improving mid-frequency contrast in sparse aperture optical imaging systems based upon the Golay-9 array,” Opt. Express 18(5), 4417–4427 (2010). [CrossRef] [PubMed]

6. N. J. Miller, M. P. Dierking, and B. D. Duncan, “Optical sparse aperture imaging,” Appl. Opt. 46(23), 5933–5943 (2007). [CrossRef] [PubMed]

7. L. M. Mugnier, G. Rousset, and F. Cassaing, “Aperture configuration optimality criterion for phased arrays of optical telescopes,” J. Opt. Soc. Am. A 13(12), 2367–2374 (1996). [CrossRef]

8. I. Tcherniavski and M. Kahrizi, “Optimization of the optical sparse array configuration,” Opt. Eng. 44(10), 103201 (2005). [CrossRef]

9. P. S. Salvaggio, J. R. Schott, and D. M. McKeown, “Genetic apertures: an improved sparse aperture design framework,” Appl. Opt. 55(12), 3182–3191 (2016). [CrossRef] [PubMed]

10. K. D. Bell, R. H. Boucher, R. Vacek, and M. Hopkins, “Assessment of large aperture lightweight imaging concepts,” in Aerospace Applications Conference, 1996. Proceedings (1996), pp. 187–203. [CrossRef]

11. J. R. Fienup and J. J. Miller, “Comparison of reconstruction algorithms for images from sparse-aperture systems,” Proceedings of SPIE - The International Society for Optical Engineering 4792, 1–8 (2002). [CrossRef]

12. D. Wang and S. Tao, “Experimental study on imaging and image restoration of optical sparse aperture systems,” Opt. Eng. 46(10), 13053–13057 (2007). [CrossRef]

13. Z. Zhou, D. Wang, and Y. Wang, “Effect of noise on the performance of image restoration in an optical sparse aperture system,” J. Opt. 13(7), 075502 (2011). [CrossRef]

14. E. Huggins, Introduction to Fourier Optics (McGraw-Hill, 1968).

15. P. S. Salvaggio, J. R. Schott, and D. M. McKeown, “Validation of modeled sparse aperture post-processing artifacts,” Appl. Opt. 56(4), 761–770 (2017). [CrossRef] [PubMed]

16. Y. Jiang, Z. Wang, M. Zhao, and L. Zhang, “MTF-based Research on the Minimum Fill Factor in Optical Sparse Aperture System,” in Fourth International Conference on Intelligent Computation Technology and Automation(2011), pp. 612–615.

17. H. Hotelling, “Analysis of a complex of statistical variables into principal components,” Br. J. Educ. Psychol. 24(6), 417–520 (1933). [CrossRef]

18. K. Elsayed and C. Lacor, “Comparison of Kriging, RBFNN, RBF and polynomial regression surrogates in design optimization,” in Eleventh International Conference of Fluid Dynamics(2013).

19. B. Hunt, “A matrix theory proof of the discrete convolution theorem,” IEEE Trans. Audio Electroacoust. 19(4), 285–288 (1971). [CrossRef]

20. B. Hunt, “A Theorem on the Difficulty of Numerical Deconvolution,” AC-20, 94–95 (1972). [CrossRef]

21. B. R. Hunt, “The Application of Constrained Least Squares Estimation to Image Restoration by Digital Computer,” IEEE Trans. Comput. 22(9), 805–812(1973).

22. C. R. Vogel, “Computational Methods for Inverse Problems,” Frontiers in Applied Mathematics (2002).

Image	Noise	SIV ISNR	SV ISNR
cameraman	20 dB	2.087	2.107
	30 dB	3.644	3.799
	50 dB	5.517	5.903
house	20 dB	2.797	3.055
	30 dB	4.426	5.062
	50 dB	6.382	7.747
chemical plant	20 dB	3.086	3.430
	30 dB	4.485	5.247
	50 dB	5.593	6.797

Image	Noise	SIV ISNR	SV ISNR
cameraman	20 dB	2.087	2.107
	30 dB	3.644	3.799
	50 dB	5.517	5.903
house	20 dB	2.797	3.055
	30 dB	4.426	5.062
	50 dB	6.382	7.747
chemical plant	20 dB	3.086	3.430
	30 dB	4.485	5.247
	50 dB	5.593	6.797

Mid-frequency MTF compensation of optical sparse aperture system

Abstract

1. Introduction

2. Principle of mid-frequency MTF decreasing and missing of OSA system

3. Mid-frequency MTF compensation for large filling factor OSA system

3.1. Decomposition and acquisition of Space Variant PSF

3.1.1. Measurement of PSFs

3.1.2. Unknown PSFs acquisition and storage by decomposition and interpolation

3.2. The SV image restoration algorithm

3.2.1. Calculation of deconvolution

3.2.2. Regularization of deconvolution problem

3.2.3. Solving of algebraic equation

3.3. Simulation experiment

4. Mid-frequency MTF compensation for small filling factor OSA system

4.1 Method of compensating missing mid-frequency information of OSA system

4.2 Feasibility of compensation of missing mid-frequency MTF

4.3 MTF analysis of OSA system and single aperture system

4.4 Image fusion method

4.5 Simulation experiments

5. Summary and conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (15)

Tables (2)

Equations (42)

Optics Express