Optik 124 (2013) 5542–5547
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.de/ijleo
Fourier Telescopy extrapolation based on the sampling theorem Zhisheng Zhou a,∗ , Bin Xiangli a,b , Wenxi Zhang b , Yang Li b , Xinxin Kong b a b
Department of Precision Machinery and Precision Instrumentation, University of Science and Technology of China, Hefei, Anhui 230027, China Academy of Opto-Electronics, Chinese Academy of Sciences, Beijing 100094, China
a r t i c l e
i n f o
Article history: Received 31 October 2012 Accepted 30 March 2013
Keywords: Fourier Telescopy Circular array Extrapolation Sampling theorem Resolution
a b s t r a c t Fourier Telecopy (FT) technique plays a promising role in deep space target detection. To improve the imaging quality of circular transmitting telescope array FT, a new processing method which adopts an extrapolation method based on the sampling theorem is proposed. Simulation results show that the extrapolation method can obviously improve system imaging resolution as well as reconstructed image quality. Even at a relative low SNR level the method improves system imaging resolution to a certain extent. © 2013 Elsevier GmbH. All rights reserved.
1. Introduction Fourier Telescopy (FT) is a computational imaging technique, which can create images of a distant dim object with high resolution [1–3], and it has been studied extensively for deep space target imaging under the USA Air Force programs [4–7]. Since FT is a synthetic array imaging technique, the arrangement of the transmitting telescope array is important. Different types of transmitting telescope arrays have their own advantages and disadvantages which results in different imaging quality. Cursory experience has shown that best Fourier inversion occurs for a given number of frequency samples when frequency sampling is uniform (Holmes, OSA Annual Meeting). Though sampling of some types of arrays such as circular array and rotating linear array is not uniform, these arrays are more competitive in some cases when takes engineering practical condition into account. Under the program SAINT a transmitting array with not uniform frequency sampling which looks like a circular array to some extent is proposed for more distinct frequency samples [6,7]. So to improve the imaging quality for a given not uniform array is interesting and meaningful. However related work has not been reported according to the authors’ knowledge. As a type of not uniform sampling array, circular array has the advantages of compact structure and more distinct frequency samples, but its highest measured frequency as well as data efficiency is low. Taking advantage of large numbers of distinct frequency
∗ Corresponding author at: Room 1102, Main Building of Academy of OptoElectronics, No. 9, DengZhuangNan Road, Haidian District, Beijing city, Beijing 100094, China. Tel.: +86 18210129935/01082178622; fax: +86 01082178674. E-mail address:
[email protected] (Z. Zhou). 0030-4026/$ – see front matter © 2013 Elsevier GmbH. All rights reserved. http://dx.doi.org/10.1016/j.ijleo.2013.03.148
samples, an extrapolation method based on the sampling theorem can be applied to obtaining uniform frequency samples with higher highest measured frequency and higher data efficiency from original samples. The imaging resolution as well as reconstructed image quality is obviously improved after the extrapolation. The rest of the paper is organized as follows. Section 2 describes the basic theory of the extrapolation method based on the sampling theorem. In Section 3, circular array FT extrapolation is introduced and some details are discussed. Section 4 gives out simulation results and practical limitation analyses are performed in Section 5. Section 6 concludes the paper. 2. Extrapolation based on the sampling theorem [8,9] Suppose that a two-dimensional object with intensity distribution Ig (x, y) is bounded to the region (−L/2, L/2) on the x axis as well as y axis. By the Whittaker–Shannon sampling theorem, the object spectrum Gg (u, v) can be written in terms of its sample values at frequencies n/L Gg (u, v) =
∞ ∞
Gg (m/L, n/L) sin c[L(u − m/L)] sin c[L(v − n/L)]
m=−∞n=−∞
(1) Due to the limited passband of the optical system, values of Gg (m/L, n/L) can be found for only a few low-integer values of m, n. Suppose −N ≤ m ≤ N, −N ≤ n ≤ N, so that the approximation Gg (u, v) ≈
N N
Gg (m/L, n/L) sin c[L(u−m/L)] sin c[L(v−n/L)]
m=−N n=−N
(2)
Z. Zhou et al. / Optik 124 (2013) 5542–5547
Fig. 1. Schematic of circular array transmitting telescope array FT.
5543
Fig. 2. Frequency distributions of circular array and uniform array FT (for interpretation of the references to color in text, the reader is referred to the web version of the article).
To determine the sample values outside the observable passband, measure the values of Gg (u, v) at any (2N + 1)2 distinct frequencies (ui , vi ) within the passband. Thus it generates a set of (2N + 1)2 simultaneous equations
3. Circular array FT extrapolation
star) is also shown for comparison. To ensure the validity of the comparison, both arrays have 19 transmitting telescopes and the same separation distance between adjacent telescopes. Circular array can achieve 342 distinct frequency samples while uniform sampling array can only achieve 169 distinct samples. Nevertheless, the highest measured frequency of circular array is lower than that of uniform sampling array (in 45◦ direction), and consequently its resolution is lower. Furthermore, the gap is more and more obvious as transmitting telescope increases. Another problem of circular array is that measured frequencies tend to be more and more concentrated as they get higher, which results in the fact that lots of high frequencies crowd in a narrow annular space and thus measured spectra data are of low efficiency. Because of the above weaknesses, circular array has inferior imaging quality comparing with uniform sampling array though it obtains more measured spectra and computer simulations have proven this inference. The above-mentioned extrapolation method is suitable for circular array FT. First, the extrapolation method deals with data in frequency domain while FT obtains spectra of target image through demodulation, there is no need to make data domain transform. Second, the more distinct measured samples, the higher frequency beyond the passband can be achieved through the extrapolation, and meanwhile circular array FT obtains large numbers of distinct frequency samples. Third, the extrapolation obtains efficient uniform sampling spectra from low-efficiency circular array spectra. So it is theoretically a good way to improve circular array FT imaging quality. It is need to be mentioned that though this paper focuses on circular array, any not uniform transmitting telescope array with the same characteristics of measured values distribution can apply the extrapolation method.
3.1. Circular array FT
3.2. Resolution improvement
A schematic of circular transmitting telescope array FT is shown in Fig. 1. N transmitting telescopes are arranged to a regular polygon and each vector connecting any pair of transmitting telescope positions is proportional to the corresponding measured spatial frequency of target image [6]. When N is an odd number, any two vectors connecting telescope positions are not the same, thus all measured frequencies are distinct. Circular array has the advantages of compact structure and more distinct frequency samples for a given number transmitting telescopes comparing with uniform sampling array, however, its highest measured frequency and data efficiency is lower. Fig. 2 shows measured frequency distribution of circular array (red point), and frequency distribution of uniform sampling array (blue
Circular array FT can obtain large numbers of distinct frequency samples, which means N can be very large and the highest frequency N/L figured out through the extrapolation is likely to exceed original frequencies, as a result the system imaging resolution is likely to improve. Assume that the number of circular array FT transmitting telescopes is M, and the separation distance between adjacent telescopes is D. It is easy to calculate the largest distance Dmax between telescopes
Gg (ui , vi ) ≈
N N
Gg (m/L, n/L) sin c[L(ui − m/L)] sin
m=−N n=−N
c[L(vi − n/L)] i = 1, 2 · · · (2N + 1)2
(3)
This is a set of (2N + 1)2 linear equations in (2N + 1)2 unknown Gg (m/L, n/L), and the set of equations can be represented by the simple matrix equation
=
g
(2N+1)2 ×1
C
g
(2N+1)2 ×(2N+1)2 (2N+1)2 ×1
(4)
Column vector g consists of the (2N + 1)2 measured values Gg (ui , vi ) and column vector g consists of the (2N + 1)2 unknown Gg (m/L, n/L). Coefficient matrix C is composed of corresponding sinc·sinc functions. It is possible to prove that as long as the measured frequencies (ui , vi ) are distinct, the determinant of C is nonzero, and therefore the inverse matrix C−1 exists. Thus in principle the unknown g can be figured out through the equation
=
g 2
(2N+1) ×1
C−1
g
(2N+1)2 ×(2N+1)2 (2N+1)2 ×1
(5)
Based on the method, the more measured frequencies, the larger N is, and thus the more spectrum information beyond the limited passband can be achieved.
Dmax = D
sin[(M − 1)/2M] sin(/M)
(6)
5544
Z. Zhou et al. / Optik 124 (2013) 5542–5547
Table 1 Fe /Fc of circular array FT of certain system parameters. (nm)
M
D (m)
R (km)
L (m)
Fe /Fc
1064
41
0.3
1000
2/3/4
2.6/1.7/1.3
The wavelength of transmitting laser beams is and the range to the target is R, so the highest measured frequency Fc on the u axis is defined by Fc =
Dmax R
(7)
The number of distinct distances between telescopes is (M − 1)/2, and for every distance there are M orientations in the upper half quadrant, so the number of measured spectrum samples is (M − 1)/2 × M × 2 including the conjugate items. Select (2N + 1)2 measured samples and solve the Eq. (4) to figure out the unknown Gg (m/L, n/L), where L is the size of target. The highest sampling frequency Fe on the u axis is defined by Fe =
Nmax (2Nmax + 1)2 ≤ (M − 1) × M L
(8)
It is generally believed that resolution is proportional to the highest frequency that a system can achieve, thus the value Fe /Fc is used to be an index of resolution estimate. Table 1 shows Fe /Fc of circular array FT of certain system parameters. Calculated result shows that for different size targets, the extrapolation method can improve resolution theoretically. However, the size of target is not known accurately in practice. In order to completely reconstruct the object spectra, the sample distance 1/Ls must obey the Whittaker–Shannon sampling theorem that 1/Ls < 1/L, as a result the highest frequency Nmax /Ls figured out is lower than Nmax /L and Fe /Fc is less than calculated result. 3.3. TSVD method In mathematics as long as the measured frequencies (ui , vi ) are distinct, the inverse matrix C−1 exists and Eq. (4) can be solved accurately. But actually C is very ill-conditioned and Eq. (4) is an ill-conditioned linear system, so it is hard to achieve inverse matrix C−1 and solution directly from Eq. (5) varies from real solution greatly. Truncated SVD (TSVD) method is a relatively effective solution to ill-conditioned linear system [10–12]. First carry out singular value decomposition (SVD) as follows 2
C
2
(2N+1) ×(2N+1)
=
Fig. 3. Distributions of frequency of original and extrapolation (for interpretation of the references to color in text, the reader is referred to the web version of the article).
For original circular array measured values, reconstruct formula is defined as follow
(M−1)/2 2M
i(x, y) =
Gmn exp
m=1
Gmn = G
j
2 Dm [cos(n )x + sin(n )y] R
n (Dm 2 − Dm−1 2 )/2
n=1
Dm Dm cos(n ), sin(n ) R R
(11) where Gnm is the measured spectrum at the frequency proportional to the vector with length Dm and angle n , and n (Dm 2 − Dm−1 2 )/2 is the division area for integral in spectrum domain. Actually this formula is a Riemann integral of Fourier inverse transform. For values figured out through extrapolation, reconstruct formula is defined as follow i(x, y) =
N N
2
Gg (m/L, n/L) exp j
m=−N n=−N
R
(mx/L + ny/L)
(12)
It is also a Riemann integral of Fourier inverse transform, but there is no need to consider the division area since for uniform sampling values the division area is the same. 4. Simulation results
U
2
2
D
2
2
2
V
2
(9)
(2N+1) ×(2N+1) (2N+1) ×(2N+1) (2N+1) ×(2N+1)
4.1. Targets used in simulations
where U and V is unitary matrix, D is the singular value matrix in which singular values are arranged in descending order. TSVD method abandons the unreliable small singular values and keeps the reliable large ones. Truncated parameter is based on many kinds of rules. The kept singular values form a new matrix Dn and then approximate solution g¯ of unknown g is defined by T g¯ = VD−1 n U g
Computer simulations are carried out to test the extrapolation method, and Fig. 4 shows the images of targets used in simulations. The target image for resolution estimate is composed by three lines with equal width and equal separation distance (a), while the target
(10)
3.4. Reconstruction By the extrapolation method uniform sampling values at frequencies n/L can be figured out from original circular array measured values. Fig. 3 shows distributions of original measured values (blue point) as well as values after extrapolation (red star) with system parameters described in Table 1, of which the target size is 3 m. Then to conduct Fourier inverse transform to reconstruct target image.
Fig. 4. Target images used in simulations.
Z. Zhou et al. / Optik 124 (2013) 5542–5547
image for reconstructed image quality estimate is a grayscale image of a static satellite which is often used as target in FT simulations.
5545
extrapolation, satellite contours and small details are much clearer in reconstructed image especially for the 4 m target.
4.2. Resolution improvement First a simulation is conducted to verify resolution improvement. The system parameters are as described in Table 1 and in principle the highest measured frequency Fc is 3.68/m. It is generally believed that the minimum separation distance between line centers that the system can distinguish is 1/Fc (0.27 m), thus the separation distance is chosen to be 0.27 m first. Fig. 5(a) is the reconstructed image of measured circular array spectra and Fig. 5(b) is the reconstructed image after extrapolation. It is shown that as predicted system resolution is about 0.27 m and the lines are recognized in both directly reconstructed image and extrapolation image, however the lines in the latter are more clear-cut. Fig. 5(c) and (d) shows result of 0.15 m separation distance. Since 0.15 m is far beyond the system resolution, it is not possible to distinguish the three lines in directly reconstructed image, but as Fig. 5(d) shown it is still easy to distinguish three lines after extrapolation. Fig. 5(e) and (f) shows result of 0.10 m separation distance. Though the three lines can be recognized after extrapolation, the edges are so blurry that it is reasonable to conclude that 0.1 m is nearly the resolution the method can achieve. Anyway there is no doubt that the extrapolation method can improve resolution obviously according to the simulation results. 4.3. Reconstructed image quality Reconstructed image quality is another important evaluation criterion of an imaging system. Here Image Strehl ratio is used as the image quality metric [13]. The system parameters are also as described in Table 1. Fig. 6(a), (c) and (e) shows images directly reconstructed from measured circular array spectra and the target size is 2 m, 3 m and 4 m, respectively. For the 2 m target measured high frequencies are not enough high so that small details of reconstructed satellite are serious missing; on the contrary, for the 4 m target measured low frequencies are not enough low so that reconstructed satellite contours are not clear. The Image Strehl ratio of reconstructed image and target image is 0.86, 0.92 and 0.82, respectively. Fig. 6(b), (d) and (f) shows corresponding reconstructed images after extrapolation and the image quality visibly has been improved. The Image Strehl ratio is advanced to 0.93, 0.94 and 0.96 accordingly. Though Image Strehl ratio does not vary much after
5. Practical limitation analyses Simulations in Section 4 do not consider noise, but all methods for extrapolating bandwidth beyond the diffraction limit are known to be extremely sensitive to noise in the measured data [8,14,15]. As it is mentioned in Section 3.3, matrix C is very ill-conditioned. Moreover, as attempts are made to estimate more and more values of the spectra on the sampling points outside the observable passband, the matrix becomes more and more ill-conditioned. The growth of noise sensitivity is extremely rapid. Thus to analyze FT extrapolation effect in the case of noise is necessary and essential. In order to test the anti-noise performance of FT extrapolation, Gaussian white noise is added to the measured spectra demodulated from echo signals. Define signal-noise ratio (SNR) as the ratio of spectrum signal power to Gaussian white noise power. In reality the SNR level in laboratory circumstance is about 50–500 while in field circumstance is about 50–150. Simulations are conducted to see resolution performance at different SNR level. Fig. 7(a) and (b) shows simulation result of the three lines target image with 0.27 m separation distance and 50 SNR level. Although extrapolation result (b) is certainly not as good as Fig. 5(b) which is an ideal extrapolation result, it still displays clearer three lines than directly reconstructed result (a) does. Fig. 7(c) and (d) shows simulation result with 0.23 m separation distance and 50 SNR level. Obviously 0.23 m is beyond system distinguish ability (c) but it is still easy to recognize three lines after extrapolation (d). However as the separation distance decreases, reconstructed result gets worse rapidly (e and f), thus it is reasonable to conclude that 0.23 m is the resolution of extrapolation method at 50 SNR level. Fig. 8 shows resolution change of extrapolation method at different SNR level (red star). Comparing with resolution of directly reconstructed method (blue point), resolution after extrapolation method is improved visibly, but it is far lower than resolution after extrapolation method without noise (magenta point). It means on the one hand the extrapolation method has certain anti-noise performance, even in the situation of low SNR level it improves imaging resolution; on the other hand its anti-noise performance is not very strong, low SNR level badly restrains its ability to improve imaging resolution.
Fig. 5. (a) 0.27 m directly, (b) 0.27 m extrapolation, (c) 0.15 m directly, (d) 0.15 m extrapolation, (e) 0.10 m directly and (f) 0.10 m extrapolation.
5546
Z. Zhou et al. / Optik 124 (2013) 5542–5547
Fig. 6. (a) 2 m directly, (b) 2 m extrapolation, (c) 3 m directly, (d) 3 m extrapolation, (e) 4 m directly and (f) 4 m extrapolation.
Fig. 7. (a) 0.27 m directly, (b) 0.27 m extrapolation, (c) 0.23 m directly, (d) 0.23 m extrapolation, (e) 0.20 m directly and (f) 0.20 m extrapolation.
To improve the anti-noise performance, efforts should be made to try matrix pretreatment to refine the morbidity of matrix C, and other algorithm for solving ill-conditioned linear system should be test. All these work remains to be done in future.
6. Conclusions
Fig. 8. Resolution versus different SNR level (for interpretation of the references to color in text, the reader is referred to the web version of the article).
To improve the imaging quality of circular transmitting telescope array FT, a new processing method which adopts an extrapolation method based on the sampling theorem is proposed. The FT extrapolation method makes full use of distribution characteristics of measured spectra of circular array and obtains high-efficiency uniform sampling spectra with higher highest measured frequency, which results in the facts that it improves circular array FT imaging resolution and reconstructed image quality theoretically. Computer simulation results in fact strongly support the conclusion. Moreover, simulation results show that even in very low SNR level the extrapolation method improves system imaging resolution. Therefore, the proposed processing method is meant to
Z. Zhou et al. / Optik 124 (2013) 5542–5547
be a good way to improve imaging quality of not uniform sampling transmitting telescope array FT. Acknowledgements This work was supported by the National High Technology Research and Development Program of China. The authors would like to thank the program executive committee for the continuous supports. References [1] P.A. Bakut, Theoretical Studies of Fourier Telescopy for Deep Space Imaging, International Information Academy, Moscow, 1999, pp. 1–44. [2] R.B. Holmes, T.J. Brinkley, Reconstruction of Images of Deep Space Objects Using Fourier Telescopy, Digital Image Recovery and Synthesis IV, vol. 3815, SPIE, Denver, Colorado, 1999, pp. 11–22. [3] I. Valery, Mandrosov, Compact schematic of high resolution Fourier Telescopy imaging in strongly inhomogeneous atmosphere, Automatic Target Recognition XII, SPIE 4726 (2002) 295–303. [4] S.D. Ford, D.G. Voelz, V.L. Gamiz, S.L. Storm, S.R. Czyzak, Geo light imaging national testbed (GLINT) past, present, and future, Digital Image Recovery and Synthesis IV, vol. 3815, SPIE, Denver, Colorado, 1999, pp. 2–10.
5547
[5] V.L. Gamiz, R.B. Holmes, S.R. Czyzak, D.G. Voelz, GLINT: program overview and potential science objectives, Imaging Technology and Telescopes, SPIE 4091 (2000) 304–315. [6] B. Spivey, J. Stapp, D. Sandler, Phase closure and object reconstruction algorithm for Fourier Telescopy applied to fast-moving targets, Unconventional Imaging II, SPIE 6307 (2) (2006) 1–16. [7] J. Stapp, B. Spivey, L. Chen, et al., Simulation of a Fourier Telescopy imaging system for objects in low earth orbit, Unconventional Imaging II, SPIE 6307 (1) (2006) 1–11. [8] J.W. Goodman, Introduction to Fourier Optics, second ed., MaGraw-Hill Companies Inc., New York, 1996, pp. 162–165. [9] D. Verhoeven, Limited-data computed tomography algorithms for the physical sciences, Appl. Opt. 32 (20) (1993) 3736–3754. [10] P.C. Hansen, S. Christiansen, An SVD analysis of linear algebraic equations derived from first kind integral equations, J. Comput. Appl. Math. 12/13 (1985) 341–357. [11] P.C. Hansen, The truncated SVD as a method for regularization, Bit Numer. Math. 27 (4) (1987) 534–553. [12] Peiliang Xu, Truncated SVD methods for discrete linear ill-posed problems, Geophys. J. Int. 135 (2) (1998) 505–514. [13] E. Louis Cuellar, J. Stapp, J. Cooper, Laboratory and field experimental demonstration of a Fourier Telescopy imaging system, Unconventional Imaging, SPIE 5896 (D) (2005) 1–15. [14] G. Toraldo Di Francia, Degrees of freedom of an image, J. Opt. Soc. Am. 59 (7) (1969) 799–804. [15] R.W. Gerchberg, Super-resolution through error energy reduction, Opt. Acta 21 (9) (1974) 709–720.