- Email: [email protected]
A compressed sensing approach for enhancing infrared imaging resolution
Optics & Laser Technology 44 (2012) 2354–2360
Contents lists available at SciVerse ScienceDirect
Optics & Laser Technology journal homepage: www.elsevier.com/locate/optlastec
A compressed sensing approach for enhancing infrared imaging resolution Long-long Xiao a,n, Kun Liu a, Da-peng Han a, Ji-ying Liu b a b
Institution of Space Technology, College of Aerospace and Material Engineering, National University of Defense Technology, Deya Road, Changsha, China College of Science, National University of Defense Technology, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 18 March 2012 Received in revised form 12 April 2012 Accepted 18 April 2012 Available online 9 May 2012
This paper presents a novel approach for improving infrared imaging resolution by the use of Compressed Sensing (CS). Instead of sensing raw pixel data, the image sensor measures the compressed samples of the observed image through a coded aperture mask placed on the focal plane of the optical system, and then the image reconstruction can be conducted from these samples using an optimal algorithm. The resolution is determined by the size of the coded aperture mask other than that of the focal plane array (FPA). The attainable quality of the reconstructed image strongly depends on the choice of the coded aperture mode. Based on the framework of CS, we carefully design an optimum mask pattern and use a multiplexing scheme to achieve multiple samples. The gradient projection for sparse reconstruction (GPSR) algorithm is employed to recover the image. The mask radiation effect is discussed by theoretical analyses and numerical simulations. Experimental results are presented to show that the proposed method enhances infrared imaging resolution significantly and ensures imaging quality. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Compressed sensing Infrared imaging resolution Coded aperture mask
1. Introduction Infrared cameras have all-weather imaging capability, which are widely used in space for ground observation, target reconnaissance and strategic surveillance. The most utilized infrared sensors are uncooled micro-bolometers due to their highly reliable, small, lightweight and low input power [1]. As far as know, the existing uncooled infrared sensors with class S specification have small array size (usually 320 240), which can not satisfy resolution requirements of special missions. Demands for large format and high resolution FPAs have historically meant even smaller pixels grown on even larger substrates. For typical materials used in infrared imagers (e.g., HgCdTe, InSb), the manufacturing process which should make material geometry at the pixel-level faces extreme challenges. Owing to infrared sensor’s processing technology, traditional methods for example diminishing the pixel-pitch are hardly to develop larger format arrays and improve imaging resolution. The ability to improve the resolution of existing FPA technology without physical reduction in pixel size, therefore, has some potentially powerful implications for the infrared imaging community. Compressed sensing (CS), which captures and represents compressible signals at a sampling rate significantly below the Nyquist rate, serves as a new framework for signal sampling and reconstruction based on signal sparsity or compressibility [2]. The
n Corresponding author. Tel.: þ86 731 84574194; þ 86 13548682547; fax: þ 86 731 84512301. E-mail address: [email protected] (L.-l. Xiao).
0030-3992/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.optlastec.2012.04.022
sampling rate depend on signal sparsity other than signal bandwidth. The projection data are considerably less voluminous than those acquired by traditional sampling methods. Many researchers have studied the application of CS in extensive fields, such as analog-to-information conversion [3], synthetic aperture radar processing [4], remote sensing imagery [5], magnetic resonance imaging [6], and astronomical data analysis [7]. The most noteworthy application is the single-pixel camera designed by Baraniuk’s team at Rice University [8]. A digital micro-mirror device (DMD) is used to form a measurement matrix. This device is a grid of micro-mirrors that reflect a proportion of incoming light beams toward the sensor (a single photodiode). Each measurement made by the photodiode is a random combination of many pixels. By taking Klog(N/K) snapshots with a different random selection of pixels each time, the single-pixel camera acquires a recognizable image with a resolution comparable to N pixels. Inspired by the idea of the signal-pixel camera, a millimeter wave high resolution imaging system based on CS was designed in Ref. [9]. In the system, the energy radiated by the observed object passes though a space filter called spatial modulation. Then it is received by an antenna. In order to reconstruct the object, several samples are needed and the number of samples is correlated with the measurement matrix size. The simulation results show that 200 samples are required for recovering a 16 16 pixels image with a good imaging quality. This sampling method makes a low imaging rate in despite of its low cost. Ref. [10] shows that coded aperture designs based on the principles of CS lead to very simple and robust camera architectures which are particularly effective for low light video settings.
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
The simulation demonstrates that coded apertures are optimal when making compressive measurements and allowing for nonlinear reconstruction. The method achieves higher spatial resolution than a pinhole and a traditional coded aperture method called Modified Uniformly Redundant Arrays (MURAs). Ref. [11] presents a tutorial to discusses compressed sensing in the context of optical imaging devices, emphasizing the practical hurdles related to building physical devices and offering suggestions for overcoming these hurdles. In the infrared camera example, a first phase modulator makes use of Fourier optics to convolve the electromagnetic field associated with the image and a random phase pattern; A second phase modulator, located at the image plane, is necessary in order for the resulting sensing matrix satisfies the RIP (restricted isometry property). The advantage of the architecture is that it modulates phase only, thus light is not blocked as it passes through coded aperture architectures. However, this camera architecture requires two SLMs (spatial light modulation). Furthermore, the sampling process refers to both Fourier transform and inverse Fourier transform, and image detector should measure the complex image field are far less common than those that measure image intensity. Inspired by the above applications, we consider whether high resolution infrared imaging can be implemented based on coded aperture masks. Fortunately, the pixel-pitch of infrared sensors is in general much larger than that of their visible-light counterparts. So it is realizable to separate each pixel to several virtual ones by a coded aperture mask. In this paper, we only place a coded aperture on the focal plane of the optical system instead of two SLMs in Ref. [11]. This coding mode is clearly a simple design and can be manufactured easily. In this model, the sensor merely needs to be capable of measuring the image intensity (as opposed to both real and imaginary parts of the complex image). Furthermore, unlike the several measurement times described in Ref. [9], a multiplexing scheme is used to obtain multiple samples. Based on these above advantages, we describe an infrared imaging system based on CS. The rest of this paper is organized as follows: Section 2 surveys the basic theory of CS briefly. Section 3 discusses the focal plane coding method, mainly including the compressive sampling approach and the image reconstruction algorithm. Section 4 discusses the mask radiation effect which may affect the signal quality. Section 5 is the numerical experiment, in which the feasibility of the designed infrared imaging system is validated. The paper is concluded in Section 6.
2. CS framework and observation mode CS is a very efficient and fast growing signal recovery framework. It firstly measures the inner products produced by the
2355
signal and a certain matrix. The measurement is a small number of random linear combinations of signal values-much smaller than the number of signal samples nominally defining it. The signal is reconstructed with good accuracy using an optimization process from these inner products called projections. The CS model is described as min:x:1
s:t:
y ¼ FUf ¼ FUCUx ¼ YUx,
ð1Þ
N
where fAR is an unknown signal, which could be sparsely represented as f ¼ CUx in an orthonormal basis C. If there are only K(K 5 N) non-zero components of x, f is defined as being Ksparse. F denotes a M N (M 5 N) matrix called measurement matrix. Y is a sensing matrix compound by F and C. yAZM is an observation vector. In order to successfully reconstruct a signal with incomplete measurements, Y must satisfy a special property called the restricted isometry property (RIP) [12]; that is, for all K-sparse xARN, a constant dkA(0,1)exists so that 2
1dk r
:Yx:2 2 :x:2
r 1þ dk
ð2Þ
Intriguingly, many kinds of random matrices meet RIP with high probability. Ref. [13,14] proved that if F is a Gaussian random matrix, Y can meet RIP in all probability. Other common feasible matrices include random binary Bernoulli distribution matrix, partial Fourier matrix, local Hadamard matrix, Toeplitz matrix and so on. There are many optimal methods to solve the problem described in Eq. (1), such as Basis Pursuit (BP) [15], Matching Pursuit (MP) [16], iterative thresholding [17], and Total Variation (TV) [18]specially for image processing. From the above model, CS contains three critical parts: signal sparsity, random observation and recovery algorithm. The emphasis of this paper is the design of a measurement matrix for focal plane aperture coding. The observation model of this system is established as: y ¼ FUf þn, where f is the image of interest, F is a measurement matrix projecting the scene onto an M-dimension set, y is an M-dimension observation vector. n is white noise associated with the observation process. The measurement matrix F can correspond to a wide variety of imaging system models. The main problem is the estimation of F in order to recover f from y, when the dimensionality of y is now significantly less than the dimensionality of f. This paper addresses the accurate reconstruction of a high resolution image f from a single low resolution observation y via compressive coded aperture imaging. In particular, we describe the design of coded aperture patterns for incoherent imaging systems which significantly improve the accuracy with which we can reconstruct f, and theoretically characterize the performance of the system within the framework of compressive sensing.
Fig. 1. Focal plane coding mask (a) pixel coding mode (b) one pixel sparse sampling.
2356
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
3. Focal plane coding for high resolution infrared imaging In general, the widely used infrared sensors have the format of 320 240 pixels, which cannot satisfy the resolution requirements of many science experiments. For the above restriction, this section develops a CS method to improve imaging resolution by designing a suitable coded aperture mask to represent the measurement matrix F as described in Section 2.
there is an m m mask sub-array of each pixel, with each subregion filled with n n sensor pixels and the total number of lenses being l l. Then we can obtain M¼(l l) (n n) measurements at a single snapshot, while the dimension of the recovered image is N¼ (m m) (n n). According to CS theory, the object can be reconstructed accurately even when MoN. The image resolution is enhanced by the ratio m/l 1. 3.2. Implementation of the sensing step
3.1. Focal plane coding and sparse sampling Fig. 1(a) gives a 12 12 coded aperture mask array, corresponding to a 4 4 pixel sensor with a 3 3 mask sub-array for each pixel. Usually a random 0–1 Bernoulli matrix is used to design the mask because the ‘0’ and ‘1’ elements can be realized by a proper physical mask design which blocks or passes the light respectively [19]. As depicted in Fig. 1(b), the white area enable the light to pass through the mask while the black area disable it. The number of ‘‘white’’ pixels is almost the same as that of ‘‘black’’ pixels. Each sensor pixel records the light intensity of the incoming light beam encoded by a 3 3 mask. In fact, these records are the compression samples of the observed object. The object can be recovered from them by an optimization algorithm. It is the focal plane array size other than the sensor pixel-pitch that determines the imaging resolution. As the former is smaller than the latter, the imaging resolution can be significantly enhanced. The measurement of each pixel can be described as Z y ¼ f ðrÞdr ð3Þ A
where y is the pixel measurement, f(r) is the light intensity at r on the focal plane, A is the pixel area. The measurement as shown in Fig. 1(b) can be encoded as Z y¼ FðrÞf ðrÞdr ð4Þ where F(r) represents coded aperture mask array of the pixel. Discretizating r, Eq.(4) can be rewritten as y ¼ FUf
In order to obtain the compressed measurement, the key point is to design a suitable measurement matrix. As discussed in Section 2, when its elements are assigned ‘0’ or ‘1’, subjecting to a uniform probability Bernoulli distribution, the measurement matrix can satisfy RIP. The matrix representation of the multiplexing coding is 2 3 2 3 j1,1 j1,2 j1,mm y1,1 y1,2 y1,nn 6 y 7 6 j2,2 j2,mm 7 y2,2 y2,nn 7 6 j2,1 6 2,1 7 6 7¼6 7 6 ^ 7 6 ^ 7 ^ ^ ^ ^ ^ ^ 4 5 4 5 jll,1 jll,2 jll,mm yll,1 yll,2 yll,nn 2 3 f 1,1 f 1,2 f 1,nn 6 f f 2,2 f 2,nn 7 6 2,1 7 7 6 6 ^ 7 ^ ^ ^ 4 5 f mm,1 f mm,2 f mm,nn ð6Þ where ji,: is a row vector of F reshaped from the coded aperture mask sub-array of the pixel corresponding to the ith lens, f:,j is a column vector of f rearranged from the image observed by the jth pixel of the corresponding sensor, and yi,:is a row vector of y reshaped from the real measurement of the sensor corresponding to the ith lens. Multiplexing is illustrated as Fig. 3 In Fig. 3, y is the output of sensor, F is the measurement matrix from rearrangement of all coded aperture mask subarrays, and f represents the observed object. Therefore, the observation mode is y¼ FUf, where f can be recovered from y by a smart reconstruction algorithm designed in the following text.
ð5Þ
where F is the row vector reshaped from F(r), f is the column vector rearranged from f(r), and y is a compressed sample of f. It requires enough samples to reconstruct f. A feasible way to obtain enough samples is multiplexing, as shown in Fig. 2. There are several lenses forming a lens array, and all the lenses’ focal planes are on the same plane. A common sensor, which is segmented to several regions (the number of regions is the same as that of lenses), is placed on the focal plane. Each lens corresponds to one sub-region of the sensor, and observes the same object. The number of lenses is the measurement times. Considering the different coding mode for each measurement, the coded aperture mode of each pixel on the same sub-region is the identical, but varies among these sub-regions. Supposed that
Fig. 2. Focal plane coding based on multiplexing.
3.3. Sparse reconstruction algorithm There are many reconstruction algorithms for recovering sparse signals. BP is the most representative algorithm, but it is not suitable for large-scale two-dimensional image problem because of its high computational complexity. Recently, greedy algorithms such as Orthogonal Matching Pursuit (OMP) and its developments provide optimal performance, guaranteeing less computational cost and storage. To make a compromise between the ability to obtain an accurate solution and retaining efficiency, the gradient projection for sparse reconstruction (GPSR) algorithm is employed in this paper [20]. Firstly, the constrained
Fig. 3. Optical multiplexing and compressed sampling base on focal plane coding.
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
2357
optimization problem as Eq. (1) can be transformed into the convex unconstrained optimization problem 2
min 12 :yYUx:2 þ l:x:1
ð7Þ
x
By splitting the variable x into its positive and negative parts ~ ¼[Y Y]. x ¼u v, and setting 1n ¼[1,1,y,1]T, z¼ [u; v] and Y Eq. (6) can be rewritten to ~ z:2 þ lU12n Uz FðzÞ min :yY 2 x
s:t:
z40
ð8Þ
The two-step gradient projection method defines its iterates zk þ 1 from the previous iterate zk as wk ¼ ðzk ak UrFðzk ÞÞ þ k
zk þ 1 ¼ zk þ b ðwk zk Þ
ð9Þ
where wk is a temporary variable, ( ) þ denotes the positivepart operator defined as (x) þ ¼max{0, x}, ak and bk are the scalar parameters confined to ak 4 0 and 0 o bk o 1. The choices of ak and bk are discussed in detail in Ref. [20]. x can be obtained from the convergence of zk þ 1. The method iterates along the negative gradient direction so that it is very fast and efficient.
4. Mask radiation effect analysis As discussed in the above section, a coded aperture mask is placed above the sensor. The mask blocks about half of the incoming light and brings diffractive effect. All the unwanted irradiance reduction which would affect the signal quality will be theoretically analyzed in this section. In Fig. 1(b), we illustrate a 3 3 mask sub-array for each pixel. Half of the number of the mask sub-array is not an integer (9/ 2¼4.5), which makes array elements unequal probability and results in degradation of the measurement precision. In the following text, all the analyses and experiments are based on m m mask sub-arrays where m is an even. 4.1. Theoretical analysis of diffraction Without diffractive effect, half of the incoming irradiance permeates the mask and is recorded by the sensor. The perfect expression is described as Eq. (5) in Section 3.1. However, diffraction refers to the behavior of an optical wave when its lateral extent is confined by an aperture. It accounts for the fact that light rays do not follow strictly rectilinear paths when the wave is disturbed on its boundaries. In our system, the size of coded aperture is in the degree of um which is on the order of the radiation wavelength, so the effects of diffraction become most apparent. Scalar diffraction can clearly be used for describing free-space optical (FSO) propagation and provide a reasonable approximation of the principle effects of the propagation. All the analyses in this paper assume scalar diffraction. Because of the short distance between the mask and the sensor (usually several ums), Fresnel diffraction and Fraunhofer diffraction are not suitable for this optical propagation. Rayleigh– Sommerfeld diffraction is the most accurate diffraction solution considered in this system [21]. Consider the propagation of monochromatic light from a 2D plane (source plane) indicated by the coordinate variables x P and Z (Fig. 4). At the source plane, an area defines the extent of a source or an illuminated aperture. The field distribution in the source plane is given by U1(x,Z), and the field U2(x,y) in a distant observation plane can be predicted using the first
Fig. 4. Propagation geometry for parallel source and observation planes.
Rayleigh–Sommerfeld diffraction solution [20]. ZZ z expðjkr 12 Þ U 2 ðx,yÞ ¼ dxdZ U 1 ðx, ZÞ jl P r 212
ð10Þ
where l is the optical wavelength; k is the wavenumber, which is equal to 2p/l for free space; z is the distance between the centers of the source and observation coordinate systems; and r12 is the distance between a position on the source plane and a position in the observation plane. x and Z are variables of integration, and the integral limits correspond to the area of the source S. With the source and observation positions defined on parallel planes, the distance r12 is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð11Þ r 12 ¼ z2 þ ðxxÞ2 þ ðyZÞ2 Expression (10) is, in general, a superposition integral, but with the source and observation areas defined on parallel planes, it becomes a convolution integral, which can be written as ZZ U 2 ðx,yÞ ¼ U 1 ðx, ZÞhðxx,yZÞdxdZ ð12Þ where the general form of the Rayleigh–Sommerfeld impulse response is hðx,yÞ ¼
z expðjkrÞ jl r2
and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ z2 þ x2 þ y2
ð13Þ
ð14Þ
The Fourier convolution theorem is applied to rewrite Eq. (12) as U 2 ðx,yÞ ¼ J1 fJfU 1 ðx,yÞgJfhðx,yÞgg ¼ J1 fJfU 1 ðx,yÞgHðf x ,f y Þg,
ð15Þ
where H is the Rayleigh–Sommerfeld transfer function given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð16Þ Hðf x ,f y Þ ¼ exp jkz 1ðlf x Þ2 ðlf y Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 f x þf y o1=l must be satisfied for propaqffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 gating field components. For f x þf y 41=l, H is negative expoStrictly speaking,
nent, the irradiance degrades fast when propagates. This wave is called evanescent wave which can not bring irradiance from aperture. 4.2. Simulation of diffractive effect In this section, we make a numerical simulation to evaluate the diffractive effect of the mask. The simulation condition is according to actual applications. Assume each pixel of the sensor coded by a 4 4 square aperture mask. The side length of the square aperture is w¼8 um. The mask is illuminated by a unit amplitude plane wave with l ¼5 um. The distance between the mask and the sensor is z¼20 um. Considering that the mask may be affected by the surrounding ones, the computation window is
2358
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
Fig. 5. Simulation results of diffractive effect. (a) Focal plane coded aperture mode. (b) Sensor plane irradiance.
set to a 3 3 pixel array and only the diffractive effect of the representative center pixel is evaluated, as shown in Fig. 5(a). The coding modes of each pixel are identical which subject to a binary uniform probability Bernoulli distribution. In Section 4.1, Fourier convolution theorem is applied to describe the Rayleigh–Sommerfeld expression. For numerical simulation, discrete sampling of the source field, sampling of the transfer function, and the periodic nature of the FFT would lead to a variety of artifacts in the propagation result. The main reasons are that the transfer function is not bandlimited and the signal cannot be adequately sampled. By the theoretical analysis in Ref. [21], the Nyquist frequency fx is half the sample rate and corresponds to the maximum spatial frequency that can be adequately represented given the interval Wx. When 0:26 um r Wx r 0:4 mm, fx is larger than the source effective bandwidth (a reasonable estimate BE1.25 cycles/um) and the transfer function H can be adequately sampled. All the simulation parameters are given out: l ¼5 um, w¼8 um, Wx ¼ 0:4 um, z¼20 um. The simulation result is show in Fig. 5(b), which appears reasonable with some destructive interference features and slight spreading beyond the initial aperture width. However, each pixel senses the all irradiance passing through its corresponding coded aperture mask. We only should contrast the total sensed irradiance of each pixel with that without diffraction. The diffractive effect of each coded aperture need not be considered respectively. By discretely sampling the source field (a unit amplitude plane wave), the irradiance of a 4 4 aperture mask is 3200. After propagation, the irradiance falling on the sensor is 3053, which is about 95.4% of the source irradiance. This percentage means only a small part of the irradiance is absorbed by the evanescent wave. If the percentage of each pixel is the identical, the imaging quality is immune to diffraction and only the contrast should be adjusted to a normal state. Fortunately, the simulation results of the nine pixels as show in Fig. 5(b) illustrate that the diffractive effect of each pixel are almost the same. Therefore, in our system with the given parameters, mask diffraction makes little effect on the imaging quality.
Fig. 6. Graph of PSNR-l-m. (for each option (m¼ 4/6/8) we look for the result with the smallest l that has a PSNR above 26).
5. Infrared camera experiment The utility of optical system designs based on CS theory for improving the resolution of infrared imaging systems is demonstrated via a simulation study in this section. The image for our numerical experiment is a 240 240 pixel infrared picture collected by a remote sensing camera. It is a suitable selection to evaluate imaging performance for that it has many objects with different features. Discrete Cosine Transform (DCT) is selected to form a sparse representation matrix. The projection matrix and
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
2359
Fig. 7. Images from the numerical experiment. (a) the original infrared image with 240 240 pixel. (b) the raw image by encoding sample with 30 30. (c) the recovered image when m ¼8 and l ¼ 6 with PSNR ¼26.37 (d) the error of the original image and the recovered one.
recovery algorithm have been designed in Section 3. Peak signalto-noise ratio (PSNR) is employed to evaluate the quality of recovered images, which is defined as PSNR ¼ 10 log10 PN
i¼1
1 ðxi x0i Þ2 =N
! ð17Þ
where xi is the normalized pixel value of the original image and x0i is the counterpart of the recovered image; the maximum amplitude is 1. N is the number of pixels. In our proposed method, there is an m m mask above each pixel. Consider the existing optical technology for a coded aperture mask process and the possible diffraction effect, m is preferred to be less than 10. Through compressive sampling, the actual sensor array required is not the 240 240 format, it depends on m, l and n (n¼240 in the analysis), and its format is n n (l l/m m). In this experiment, we set m ¼8, m ¼6 and m¼4 respectively. l determines the sensor format for a fixed m. The relationship of imaging quality (PSNR) and l is illustrated in Fig. 6. As shown in Fig. 6, for a fixed m, the larger l is, the better is the recovered image quality, but the number of the required sensor pixels n increases accordingly. An image quality with PSNR better than 26 dB can ensure that the image is suitable for surveillance or reconnaissance. We make a tradeoff between image quality and imaging resolution as labeled on the three curves in Fig. 6. For m¼8, m ¼6 and m ¼4, we set l ¼6, l¼5 and l¼ 3 correspondingly, the required sensor formats are 180 180, 200 200 and 180
180 with PSNR 26.37 dB, 27.76 dB and 26.09 dB, respectively. For the same format 180 180, the PSNR of the recovered image for m¼8 and l ¼6 is better than that for m¼4 and l ¼3. The main reason is that the projection matrix described in Section 3.2 with m¼8 performs better than m ¼4. Fig. 7 shows the simulation results when m ¼8 and l ¼6. Fig. 7(a) shows the original infrared image of the object recorded by a 240 240 pixel sensor. When m¼8 and l¼ 6, the required sensor has a smaller array format(180 180). Each pixel records the light field intensity encoded by an 8 8 mask, so the whole scene is sampled by a 30 30 sub-array at each measurement displayed in Fig. 7(b). Fig. 7(c) shows the recovered image recovered from 6 6 measurements, which is feasible for both visual observation and further image process. There is little difference between the original image and the recovered one as shown in Fig. 7(d), which also validates the feasible imaging system. The experiment shows that a 240 240 pixel image with PSNR¼26.37 is obtained by a 180 180 pixel sensor due to the implementation of the coded aperture mask, that is to say the resolution improve with a ratio 33% (m¼8 and l¼ 6). Owing to the simple physical implementation of mask, the ratio 33% is considerable.
6. Conclusion and discussion This paper describes a focal plane coded aperture mask method for improving infrared imaging resolution based on
2360
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
CS. The infrared rays of the observed object spread across the focal plane coded aperture mask and are sparsely recorded on the sensor. Optical multiplexing makes it possible to sample the observed object several times at one exposure. This mechanism guarantees the imaging quality if enough samples are properly obtained. For this system, diffractive effect makes little impact on observation. Numerical experiments show that image resolution can be enhanced by 33% without image quality degradation (PSNR 4 26.00). The advantages of this proposed method are the followings. First, the sampling strategy is simpler than others and the multiplexing enables only one exposure to be enough for recovery original images. Secondly, the light intensity is coded in spatial domain; therefore, the sensor merely needs to be capable of measuring the image intensity as opposed to both real and imaginary parts of the complex image. Thirdly, DCT is adapted for image sparse representation, which makes the recovery algorithm work in real domain without any complex computation. All these make the implementation of the proposed system easily. However, the main shortage of this sampling strategy is that the coded aperture mask is fixed which can not change in time and would not be usable in many applications. Moreover, a coded aperture mask with a zero-mean is not physically realizable in optical systems. This is incompatible with the fact that generative models for random projection matrices used in CS often involve drawing elements independently from a zero-mean probability distribution. Therefore, it is necessary to accurately model practical optical systems with some compensation, for example shifting elements, which maybe negatively impacts the performance of the proposed reconstruction algorithm. Further works include the optimization of measurement matrices and recovery algorithms with emphasis on measurement robustness and computational complexity. In fact, if each lens matches a single sensor and all lenses observe the same object, the resolution is merely related with coded aperture mask mode, and the recovery precision only depends on the number of lenses. This method can enhance the resolution by several times at the cost of increasing system complexity, which is under study in our group.
Acknowledgment This work was supported by the National 863 project of China (Grant No. 2011AAnnn053).
References [1] Murphy R, Kohin M, Backer B. Recent Developments in Uncooled IR technology, in SPIE 2000:12–6. [2] Donoho DL. Compressed Sensing,. IEEE Transactions on Information Theory 2006;52(4):1289–306. [3] Laska, JN, Kirolos, S, Duarte, MF, Ragheb, TS, Baraniuk, RG and Massoud, Y. Theory and implementation of an analog-to-information converter using random demodulation, in Proc. of IEEE International Symposium on Circuits and Systems, 2007, 1959-1962. [4] Baraniuk, R and Philippe, S. Compressive radar imaging, in Proc. of IEEE Radar Conference, 2007, 128-133. [5] Ma JW. Single-pixel remote sensing,. IEEE Geoscience and Remote Sensing Letters 2009;6(2):199–203. [6] Lustig M, Donoho DL, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magnetic Resonance In Medicine 2007;58(6): 1182–95. [7] Bobin J, Starck JL, Ottensamer R. Compressed sensing in astronomy,. IEEE Journal of Selected Topics in Signal Processing 2009;28(2):609–18. [8] Duarte M, Davenport M, Takhar D, Laska J, Sun T, Kelly K, et al. Single-pixel imaging via compressive sampling,. IEEE Signal Processing Magazine 2008;25(2):83–91. [9] Han B, Zhou CL. High-resolution infrared or millimeter wave imaging based on compressed sensing. Infrared 2010;31(8):1–5 in Chinese. [10] Marcia Roummel F, Harmany Zachary T, Willett Rebecca M. Compressive coded aperture imaging. In: Proceeding of SPIE; 2009. p. 72460G. [11] Willett Rebecca M, Marcia Roummel F, Nichols Jonathan M. Compressed sensing for practical optical imaging systems: a tutorial. Optical Engineering 2011;50(7):072601. [12] Candes E. The restricted isometry property and its implications for compressed sensing. Compte Rendus de l’Academie des Sciences 2006;346(I): 589–92. [13] Candes E, Tao T. Decoding by linear programming,. IEEE Transaction on Information Theory 2005;51(12):4203–15. [14] Candes E, Romberg J, Tao T. Stable signal recovery from incomplete and inaccurate measurements. Communications on Pure and Applied Mathematics 2006;59(8):1207–23. [15] Donoho DL. For most large underdetermined systems of linear equations, the minimal L1 norm solution is also the sparsest solution. Communications On Pure And Applied Mathematics 2006;59(6):797–829. [16] Tropp J, Gilbert A. Signal recovery from random measurements via orthogonal matching pursuit. Transactions on Information Theory 2007;53(12):4655–66. [17] Herrity, K, Gilbert, A and Tropp, JA. Sparse approximation via iterative thresholding, in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, Washington DC, USA, 2006, 624-627. [18] Candes E, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,. IEEE Transactions on Information Theory 2006;52(2):489–509. [19] Portnoy AD, Pitsianis NP, Brady DJ, Guo J, Fiddy MA, Feldman MR, et al. Thin digital imaging system using focal plane coding. Proc. of SPIE-IS&T Electronic Imaging, SPIE 2006;6065:108–15. [20] Figueiredo MA, Nowak R, Wright SJ. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems, IEEE. Journal of Selected Topics in Signal Processing 2007;1(4):586–98. [21] Voelz David. Computational Fourier Optics_A Matlab Tutorial. Washington: SPIE Press; 2010.
Contents lists available at SciVerse ScienceDirect
Optics & Laser Technology journal homepage: www.elsevier.com/locate/optlastec
A compressed sensing approach for enhancing infrared imaging resolution Long-long Xiao a,n, Kun Liu a, Da-peng Han a, Ji-ying Liu b a b
Institution of Space Technology, College of Aerospace and Material Engineering, National University of Defense Technology, Deya Road, Changsha, China College of Science, National University of Defense Technology, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 18 March 2012 Received in revised form 12 April 2012 Accepted 18 April 2012 Available online 9 May 2012
This paper presents a novel approach for improving infrared imaging resolution by the use of Compressed Sensing (CS). Instead of sensing raw pixel data, the image sensor measures the compressed samples of the observed image through a coded aperture mask placed on the focal plane of the optical system, and then the image reconstruction can be conducted from these samples using an optimal algorithm. The resolution is determined by the size of the coded aperture mask other than that of the focal plane array (FPA). The attainable quality of the reconstructed image strongly depends on the choice of the coded aperture mode. Based on the framework of CS, we carefully design an optimum mask pattern and use a multiplexing scheme to achieve multiple samples. The gradient projection for sparse reconstruction (GPSR) algorithm is employed to recover the image. The mask radiation effect is discussed by theoretical analyses and numerical simulations. Experimental results are presented to show that the proposed method enhances infrared imaging resolution significantly and ensures imaging quality. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Compressed sensing Infrared imaging resolution Coded aperture mask
1. Introduction Infrared cameras have all-weather imaging capability, which are widely used in space for ground observation, target reconnaissance and strategic surveillance. The most utilized infrared sensors are uncooled micro-bolometers due to their highly reliable, small, lightweight and low input power [1]. As far as know, the existing uncooled infrared sensors with class S specification have small array size (usually 320 240), which can not satisfy resolution requirements of special missions. Demands for large format and high resolution FPAs have historically meant even smaller pixels grown on even larger substrates. For typical materials used in infrared imagers (e.g., HgCdTe, InSb), the manufacturing process which should make material geometry at the pixel-level faces extreme challenges. Owing to infrared sensor’s processing technology, traditional methods for example diminishing the pixel-pitch are hardly to develop larger format arrays and improve imaging resolution. The ability to improve the resolution of existing FPA technology without physical reduction in pixel size, therefore, has some potentially powerful implications for the infrared imaging community. Compressed sensing (CS), which captures and represents compressible signals at a sampling rate significantly below the Nyquist rate, serves as a new framework for signal sampling and reconstruction based on signal sparsity or compressibility [2]. The
n Corresponding author. Tel.: þ86 731 84574194; þ 86 13548682547; fax: þ 86 731 84512301. E-mail address: [email protected] (L.-l. Xiao).
0030-3992/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.optlastec.2012.04.022
sampling rate depend on signal sparsity other than signal bandwidth. The projection data are considerably less voluminous than those acquired by traditional sampling methods. Many researchers have studied the application of CS in extensive fields, such as analog-to-information conversion [3], synthetic aperture radar processing [4], remote sensing imagery [5], magnetic resonance imaging [6], and astronomical data analysis [7]. The most noteworthy application is the single-pixel camera designed by Baraniuk’s team at Rice University [8]. A digital micro-mirror device (DMD) is used to form a measurement matrix. This device is a grid of micro-mirrors that reflect a proportion of incoming light beams toward the sensor (a single photodiode). Each measurement made by the photodiode is a random combination of many pixels. By taking Klog(N/K) snapshots with a different random selection of pixels each time, the single-pixel camera acquires a recognizable image with a resolution comparable to N pixels. Inspired by the idea of the signal-pixel camera, a millimeter wave high resolution imaging system based on CS was designed in Ref. [9]. In the system, the energy radiated by the observed object passes though a space filter called spatial modulation. Then it is received by an antenna. In order to reconstruct the object, several samples are needed and the number of samples is correlated with the measurement matrix size. The simulation results show that 200 samples are required for recovering a 16 16 pixels image with a good imaging quality. This sampling method makes a low imaging rate in despite of its low cost. Ref. [10] shows that coded aperture designs based on the principles of CS lead to very simple and robust camera architectures which are particularly effective for low light video settings.
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
The simulation demonstrates that coded apertures are optimal when making compressive measurements and allowing for nonlinear reconstruction. The method achieves higher spatial resolution than a pinhole and a traditional coded aperture method called Modified Uniformly Redundant Arrays (MURAs). Ref. [11] presents a tutorial to discusses compressed sensing in the context of optical imaging devices, emphasizing the practical hurdles related to building physical devices and offering suggestions for overcoming these hurdles. In the infrared camera example, a first phase modulator makes use of Fourier optics to convolve the electromagnetic field associated with the image and a random phase pattern; A second phase modulator, located at the image plane, is necessary in order for the resulting sensing matrix satisfies the RIP (restricted isometry property). The advantage of the architecture is that it modulates phase only, thus light is not blocked as it passes through coded aperture architectures. However, this camera architecture requires two SLMs (spatial light modulation). Furthermore, the sampling process refers to both Fourier transform and inverse Fourier transform, and image detector should measure the complex image field are far less common than those that measure image intensity. Inspired by the above applications, we consider whether high resolution infrared imaging can be implemented based on coded aperture masks. Fortunately, the pixel-pitch of infrared sensors is in general much larger than that of their visible-light counterparts. So it is realizable to separate each pixel to several virtual ones by a coded aperture mask. In this paper, we only place a coded aperture on the focal plane of the optical system instead of two SLMs in Ref. [11]. This coding mode is clearly a simple design and can be manufactured easily. In this model, the sensor merely needs to be capable of measuring the image intensity (as opposed to both real and imaginary parts of the complex image). Furthermore, unlike the several measurement times described in Ref. [9], a multiplexing scheme is used to obtain multiple samples. Based on these above advantages, we describe an infrared imaging system based on CS. The rest of this paper is organized as follows: Section 2 surveys the basic theory of CS briefly. Section 3 discusses the focal plane coding method, mainly including the compressive sampling approach and the image reconstruction algorithm. Section 4 discusses the mask radiation effect which may affect the signal quality. Section 5 is the numerical experiment, in which the feasibility of the designed infrared imaging system is validated. The paper is concluded in Section 6.
2. CS framework and observation mode CS is a very efficient and fast growing signal recovery framework. It firstly measures the inner products produced by the
2355
signal and a certain matrix. The measurement is a small number of random linear combinations of signal values-much smaller than the number of signal samples nominally defining it. The signal is reconstructed with good accuracy using an optimization process from these inner products called projections. The CS model is described as min:x:1
s:t:
y ¼ FUf ¼ FUCUx ¼ YUx,
ð1Þ
N
where fAR is an unknown signal, which could be sparsely represented as f ¼ CUx in an orthonormal basis C. If there are only K(K 5 N) non-zero components of x, f is defined as being Ksparse. F denotes a M N (M 5 N) matrix called measurement matrix. Y is a sensing matrix compound by F and C. yAZM is an observation vector. In order to successfully reconstruct a signal with incomplete measurements, Y must satisfy a special property called the restricted isometry property (RIP) [12]; that is, for all K-sparse xARN, a constant dkA(0,1)exists so that 2
1dk r
:Yx:2 2 :x:2
r 1þ dk
ð2Þ
Intriguingly, many kinds of random matrices meet RIP with high probability. Ref. [13,14] proved that if F is a Gaussian random matrix, Y can meet RIP in all probability. Other common feasible matrices include random binary Bernoulli distribution matrix, partial Fourier matrix, local Hadamard matrix, Toeplitz matrix and so on. There are many optimal methods to solve the problem described in Eq. (1), such as Basis Pursuit (BP) [15], Matching Pursuit (MP) [16], iterative thresholding [17], and Total Variation (TV) [18]specially for image processing. From the above model, CS contains three critical parts: signal sparsity, random observation and recovery algorithm. The emphasis of this paper is the design of a measurement matrix for focal plane aperture coding. The observation model of this system is established as: y ¼ FUf þn, where f is the image of interest, F is a measurement matrix projecting the scene onto an M-dimension set, y is an M-dimension observation vector. n is white noise associated with the observation process. The measurement matrix F can correspond to a wide variety of imaging system models. The main problem is the estimation of F in order to recover f from y, when the dimensionality of y is now significantly less than the dimensionality of f. This paper addresses the accurate reconstruction of a high resolution image f from a single low resolution observation y via compressive coded aperture imaging. In particular, we describe the design of coded aperture patterns for incoherent imaging systems which significantly improve the accuracy with which we can reconstruct f, and theoretically characterize the performance of the system within the framework of compressive sensing.
Fig. 1. Focal plane coding mask (a) pixel coding mode (b) one pixel sparse sampling.
2356
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
3. Focal plane coding for high resolution infrared imaging In general, the widely used infrared sensors have the format of 320 240 pixels, which cannot satisfy the resolution requirements of many science experiments. For the above restriction, this section develops a CS method to improve imaging resolution by designing a suitable coded aperture mask to represent the measurement matrix F as described in Section 2.
there is an m m mask sub-array of each pixel, with each subregion filled with n n sensor pixels and the total number of lenses being l l. Then we can obtain M¼(l l) (n n) measurements at a single snapshot, while the dimension of the recovered image is N¼ (m m) (n n). According to CS theory, the object can be reconstructed accurately even when MoN. The image resolution is enhanced by the ratio m/l 1. 3.2. Implementation of the sensing step
3.1. Focal plane coding and sparse sampling Fig. 1(a) gives a 12 12 coded aperture mask array, corresponding to a 4 4 pixel sensor with a 3 3 mask sub-array for each pixel. Usually a random 0–1 Bernoulli matrix is used to design the mask because the ‘0’ and ‘1’ elements can be realized by a proper physical mask design which blocks or passes the light respectively [19]. As depicted in Fig. 1(b), the white area enable the light to pass through the mask while the black area disable it. The number of ‘‘white’’ pixels is almost the same as that of ‘‘black’’ pixels. Each sensor pixel records the light intensity of the incoming light beam encoded by a 3 3 mask. In fact, these records are the compression samples of the observed object. The object can be recovered from them by an optimization algorithm. It is the focal plane array size other than the sensor pixel-pitch that determines the imaging resolution. As the former is smaller than the latter, the imaging resolution can be significantly enhanced. The measurement of each pixel can be described as Z y ¼ f ðrÞdr ð3Þ A
where y is the pixel measurement, f(r) is the light intensity at r on the focal plane, A is the pixel area. The measurement as shown in Fig. 1(b) can be encoded as Z y¼ FðrÞf ðrÞdr ð4Þ where F(r) represents coded aperture mask array of the pixel. Discretizating r, Eq.(4) can be rewritten as y ¼ FUf
In order to obtain the compressed measurement, the key point is to design a suitable measurement matrix. As discussed in Section 2, when its elements are assigned ‘0’ or ‘1’, subjecting to a uniform probability Bernoulli distribution, the measurement matrix can satisfy RIP. The matrix representation of the multiplexing coding is 2 3 2 3 j1,1 j1,2 j1,mm y1,1 y1,2 y1,nn 6 y 7 6 j2,2 j2,mm 7 y2,2 y2,nn 7 6 j2,1 6 2,1 7 6 7¼6 7 6 ^ 7 6 ^ 7 ^ ^ ^ ^ ^ ^ 4 5 4 5 jll,1 jll,2 jll,mm yll,1 yll,2 yll,nn 2 3 f 1,1 f 1,2 f 1,nn 6 f f 2,2 f 2,nn 7 6 2,1 7 7 6 6 ^ 7 ^ ^ ^ 4 5 f mm,1 f mm,2 f mm,nn ð6Þ where ji,: is a row vector of F reshaped from the coded aperture mask sub-array of the pixel corresponding to the ith lens, f:,j is a column vector of f rearranged from the image observed by the jth pixel of the corresponding sensor, and yi,:is a row vector of y reshaped from the real measurement of the sensor corresponding to the ith lens. Multiplexing is illustrated as Fig. 3 In Fig. 3, y is the output of sensor, F is the measurement matrix from rearrangement of all coded aperture mask subarrays, and f represents the observed object. Therefore, the observation mode is y¼ FUf, where f can be recovered from y by a smart reconstruction algorithm designed in the following text.
ð5Þ
where F is the row vector reshaped from F(r), f is the column vector rearranged from f(r), and y is a compressed sample of f. It requires enough samples to reconstruct f. A feasible way to obtain enough samples is multiplexing, as shown in Fig. 2. There are several lenses forming a lens array, and all the lenses’ focal planes are on the same plane. A common sensor, which is segmented to several regions (the number of regions is the same as that of lenses), is placed on the focal plane. Each lens corresponds to one sub-region of the sensor, and observes the same object. The number of lenses is the measurement times. Considering the different coding mode for each measurement, the coded aperture mode of each pixel on the same sub-region is the identical, but varies among these sub-regions. Supposed that
Fig. 2. Focal plane coding based on multiplexing.
3.3. Sparse reconstruction algorithm There are many reconstruction algorithms for recovering sparse signals. BP is the most representative algorithm, but it is not suitable for large-scale two-dimensional image problem because of its high computational complexity. Recently, greedy algorithms such as Orthogonal Matching Pursuit (OMP) and its developments provide optimal performance, guaranteeing less computational cost and storage. To make a compromise between the ability to obtain an accurate solution and retaining efficiency, the gradient projection for sparse reconstruction (GPSR) algorithm is employed in this paper [20]. Firstly, the constrained
Fig. 3. Optical multiplexing and compressed sampling base on focal plane coding.
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
2357
optimization problem as Eq. (1) can be transformed into the convex unconstrained optimization problem 2
min 12 :yYUx:2 þ l:x:1
ð7Þ
x
By splitting the variable x into its positive and negative parts ~ ¼[Y Y]. x ¼u v, and setting 1n ¼[1,1,y,1]T, z¼ [u; v] and Y Eq. (6) can be rewritten to ~ z:2 þ lU12n Uz FðzÞ min :yY 2 x
s:t:
z40
ð8Þ
The two-step gradient projection method defines its iterates zk þ 1 from the previous iterate zk as wk ¼ ðzk ak UrFðzk ÞÞ þ k
zk þ 1 ¼ zk þ b ðwk zk Þ
ð9Þ
where wk is a temporary variable, ( ) þ denotes the positivepart operator defined as (x) þ ¼max{0, x}, ak and bk are the scalar parameters confined to ak 4 0 and 0 o bk o 1. The choices of ak and bk are discussed in detail in Ref. [20]. x can be obtained from the convergence of zk þ 1. The method iterates along the negative gradient direction so that it is very fast and efficient.
4. Mask radiation effect analysis As discussed in the above section, a coded aperture mask is placed above the sensor. The mask blocks about half of the incoming light and brings diffractive effect. All the unwanted irradiance reduction which would affect the signal quality will be theoretically analyzed in this section. In Fig. 1(b), we illustrate a 3 3 mask sub-array for each pixel. Half of the number of the mask sub-array is not an integer (9/ 2¼4.5), which makes array elements unequal probability and results in degradation of the measurement precision. In the following text, all the analyses and experiments are based on m m mask sub-arrays where m is an even. 4.1. Theoretical analysis of diffraction Without diffractive effect, half of the incoming irradiance permeates the mask and is recorded by the sensor. The perfect expression is described as Eq. (5) in Section 3.1. However, diffraction refers to the behavior of an optical wave when its lateral extent is confined by an aperture. It accounts for the fact that light rays do not follow strictly rectilinear paths when the wave is disturbed on its boundaries. In our system, the size of coded aperture is in the degree of um which is on the order of the radiation wavelength, so the effects of diffraction become most apparent. Scalar diffraction can clearly be used for describing free-space optical (FSO) propagation and provide a reasonable approximation of the principle effects of the propagation. All the analyses in this paper assume scalar diffraction. Because of the short distance between the mask and the sensor (usually several ums), Fresnel diffraction and Fraunhofer diffraction are not suitable for this optical propagation. Rayleigh– Sommerfeld diffraction is the most accurate diffraction solution considered in this system [21]. Consider the propagation of monochromatic light from a 2D plane (source plane) indicated by the coordinate variables x P and Z (Fig. 4). At the source plane, an area defines the extent of a source or an illuminated aperture. The field distribution in the source plane is given by U1(x,Z), and the field U2(x,y) in a distant observation plane can be predicted using the first
Fig. 4. Propagation geometry for parallel source and observation planes.
Rayleigh–Sommerfeld diffraction solution [20]. ZZ z expðjkr 12 Þ U 2 ðx,yÞ ¼ dxdZ U 1 ðx, ZÞ jl P r 212
ð10Þ
where l is the optical wavelength; k is the wavenumber, which is equal to 2p/l for free space; z is the distance between the centers of the source and observation coordinate systems; and r12 is the distance between a position on the source plane and a position in the observation plane. x and Z are variables of integration, and the integral limits correspond to the area of the source S. With the source and observation positions defined on parallel planes, the distance r12 is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð11Þ r 12 ¼ z2 þ ðxxÞ2 þ ðyZÞ2 Expression (10) is, in general, a superposition integral, but with the source and observation areas defined on parallel planes, it becomes a convolution integral, which can be written as ZZ U 2 ðx,yÞ ¼ U 1 ðx, ZÞhðxx,yZÞdxdZ ð12Þ where the general form of the Rayleigh–Sommerfeld impulse response is hðx,yÞ ¼
z expðjkrÞ jl r2
and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ z2 þ x2 þ y2
ð13Þ
ð14Þ
The Fourier convolution theorem is applied to rewrite Eq. (12) as U 2 ðx,yÞ ¼ J1 fJfU 1 ðx,yÞgJfhðx,yÞgg ¼ J1 fJfU 1 ðx,yÞgHðf x ,f y Þg,
ð15Þ
where H is the Rayleigh–Sommerfeld transfer function given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð16Þ Hðf x ,f y Þ ¼ exp jkz 1ðlf x Þ2 ðlf y Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 f x þf y o1=l must be satisfied for propaqffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 gating field components. For f x þf y 41=l, H is negative expoStrictly speaking,
nent, the irradiance degrades fast when propagates. This wave is called evanescent wave which can not bring irradiance from aperture. 4.2. Simulation of diffractive effect In this section, we make a numerical simulation to evaluate the diffractive effect of the mask. The simulation condition is according to actual applications. Assume each pixel of the sensor coded by a 4 4 square aperture mask. The side length of the square aperture is w¼8 um. The mask is illuminated by a unit amplitude plane wave with l ¼5 um. The distance between the mask and the sensor is z¼20 um. Considering that the mask may be affected by the surrounding ones, the computation window is
2358
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
Fig. 5. Simulation results of diffractive effect. (a) Focal plane coded aperture mode. (b) Sensor plane irradiance.
set to a 3 3 pixel array and only the diffractive effect of the representative center pixel is evaluated, as shown in Fig. 5(a). The coding modes of each pixel are identical which subject to a binary uniform probability Bernoulli distribution. In Section 4.1, Fourier convolution theorem is applied to describe the Rayleigh–Sommerfeld expression. For numerical simulation, discrete sampling of the source field, sampling of the transfer function, and the periodic nature of the FFT would lead to a variety of artifacts in the propagation result. The main reasons are that the transfer function is not bandlimited and the signal cannot be adequately sampled. By the theoretical analysis in Ref. [21], the Nyquist frequency fx is half the sample rate and corresponds to the maximum spatial frequency that can be adequately represented given the interval Wx. When 0:26 um r Wx r 0:4 mm, fx is larger than the source effective bandwidth (a reasonable estimate BE1.25 cycles/um) and the transfer function H can be adequately sampled. All the simulation parameters are given out: l ¼5 um, w¼8 um, Wx ¼ 0:4 um, z¼20 um. The simulation result is show in Fig. 5(b), which appears reasonable with some destructive interference features and slight spreading beyond the initial aperture width. However, each pixel senses the all irradiance passing through its corresponding coded aperture mask. We only should contrast the total sensed irradiance of each pixel with that without diffraction. The diffractive effect of each coded aperture need not be considered respectively. By discretely sampling the source field (a unit amplitude plane wave), the irradiance of a 4 4 aperture mask is 3200. After propagation, the irradiance falling on the sensor is 3053, which is about 95.4% of the source irradiance. This percentage means only a small part of the irradiance is absorbed by the evanescent wave. If the percentage of each pixel is the identical, the imaging quality is immune to diffraction and only the contrast should be adjusted to a normal state. Fortunately, the simulation results of the nine pixels as show in Fig. 5(b) illustrate that the diffractive effect of each pixel are almost the same. Therefore, in our system with the given parameters, mask diffraction makes little effect on the imaging quality.
Fig. 6. Graph of PSNR-l-m. (for each option (m¼ 4/6/8) we look for the result with the smallest l that has a PSNR above 26).
5. Infrared camera experiment The utility of optical system designs based on CS theory for improving the resolution of infrared imaging systems is demonstrated via a simulation study in this section. The image for our numerical experiment is a 240 240 pixel infrared picture collected by a remote sensing camera. It is a suitable selection to evaluate imaging performance for that it has many objects with different features. Discrete Cosine Transform (DCT) is selected to form a sparse representation matrix. The projection matrix and
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
2359
Fig. 7. Images from the numerical experiment. (a) the original infrared image with 240 240 pixel. (b) the raw image by encoding sample with 30 30. (c) the recovered image when m ¼8 and l ¼ 6 with PSNR ¼26.37 (d) the error of the original image and the recovered one.
recovery algorithm have been designed in Section 3. Peak signalto-noise ratio (PSNR) is employed to evaluate the quality of recovered images, which is defined as PSNR ¼ 10 log10 PN
i¼1
1 ðxi x0i Þ2 =N
! ð17Þ
where xi is the normalized pixel value of the original image and x0i is the counterpart of the recovered image; the maximum amplitude is 1. N is the number of pixels. In our proposed method, there is an m m mask above each pixel. Consider the existing optical technology for a coded aperture mask process and the possible diffraction effect, m is preferred to be less than 10. Through compressive sampling, the actual sensor array required is not the 240 240 format, it depends on m, l and n (n¼240 in the analysis), and its format is n n (l l/m m). In this experiment, we set m ¼8, m ¼6 and m¼4 respectively. l determines the sensor format for a fixed m. The relationship of imaging quality (PSNR) and l is illustrated in Fig. 6. As shown in Fig. 6, for a fixed m, the larger l is, the better is the recovered image quality, but the number of the required sensor pixels n increases accordingly. An image quality with PSNR better than 26 dB can ensure that the image is suitable for surveillance or reconnaissance. We make a tradeoff between image quality and imaging resolution as labeled on the three curves in Fig. 6. For m¼8, m ¼6 and m ¼4, we set l ¼6, l¼5 and l¼ 3 correspondingly, the required sensor formats are 180 180, 200 200 and 180
180 with PSNR 26.37 dB, 27.76 dB and 26.09 dB, respectively. For the same format 180 180, the PSNR of the recovered image for m¼8 and l ¼6 is better than that for m¼4 and l ¼3. The main reason is that the projection matrix described in Section 3.2 with m¼8 performs better than m ¼4. Fig. 7 shows the simulation results when m ¼8 and l ¼6. Fig. 7(a) shows the original infrared image of the object recorded by a 240 240 pixel sensor. When m¼8 and l¼ 6, the required sensor has a smaller array format(180 180). Each pixel records the light field intensity encoded by an 8 8 mask, so the whole scene is sampled by a 30 30 sub-array at each measurement displayed in Fig. 7(b). Fig. 7(c) shows the recovered image recovered from 6 6 measurements, which is feasible for both visual observation and further image process. There is little difference between the original image and the recovered one as shown in Fig. 7(d), which also validates the feasible imaging system. The experiment shows that a 240 240 pixel image with PSNR¼26.37 is obtained by a 180 180 pixel sensor due to the implementation of the coded aperture mask, that is to say the resolution improve with a ratio 33% (m¼8 and l¼ 6). Owing to the simple physical implementation of mask, the ratio 33% is considerable.
6. Conclusion and discussion This paper describes a focal plane coded aperture mask method for improving infrared imaging resolution based on
2360
L.-l. Xiao et al. / Optics & Laser Technology 44 (2012) 2354–2360
CS. The infrared rays of the observed object spread across the focal plane coded aperture mask and are sparsely recorded on the sensor. Optical multiplexing makes it possible to sample the observed object several times at one exposure. This mechanism guarantees the imaging quality if enough samples are properly obtained. For this system, diffractive effect makes little impact on observation. Numerical experiments show that image resolution can be enhanced by 33% without image quality degradation (PSNR 4 26.00). The advantages of this proposed method are the followings. First, the sampling strategy is simpler than others and the multiplexing enables only one exposure to be enough for recovery original images. Secondly, the light intensity is coded in spatial domain; therefore, the sensor merely needs to be capable of measuring the image intensity as opposed to both real and imaginary parts of the complex image. Thirdly, DCT is adapted for image sparse representation, which makes the recovery algorithm work in real domain without any complex computation. All these make the implementation of the proposed system easily. However, the main shortage of this sampling strategy is that the coded aperture mask is fixed which can not change in time and would not be usable in many applications. Moreover, a coded aperture mask with a zero-mean is not physically realizable in optical systems. This is incompatible with the fact that generative models for random projection matrices used in CS often involve drawing elements independently from a zero-mean probability distribution. Therefore, it is necessary to accurately model practical optical systems with some compensation, for example shifting elements, which maybe negatively impacts the performance of the proposed reconstruction algorithm. Further works include the optimization of measurement matrices and recovery algorithms with emphasis on measurement robustness and computational complexity. In fact, if each lens matches a single sensor and all lenses observe the same object, the resolution is merely related with coded aperture mask mode, and the recovery precision only depends on the number of lenses. This method can enhance the resolution by several times at the cost of increasing system complexity, which is under study in our group.
Acknowledgment This work was supported by the National 863 project of China (Grant No. 2011AAnnn053).
References [1] Murphy R, Kohin M, Backer B. Recent Developments in Uncooled IR technology, in SPIE 2000:12–6. [2] Donoho DL. Compressed Sensing,. IEEE Transactions on Information Theory 2006;52(4):1289–306. [3] Laska, JN, Kirolos, S, Duarte, MF, Ragheb, TS, Baraniuk, RG and Massoud, Y. Theory and implementation of an analog-to-information converter using random demodulation, in Proc. of IEEE International Symposium on Circuits and Systems, 2007, 1959-1962. [4] Baraniuk, R and Philippe, S. Compressive radar imaging, in Proc. of IEEE Radar Conference, 2007, 128-133. [5] Ma JW. Single-pixel remote sensing,. IEEE Geoscience and Remote Sensing Letters 2009;6(2):199–203. [6] Lustig M, Donoho DL, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magnetic Resonance In Medicine 2007;58(6): 1182–95. [7] Bobin J, Starck JL, Ottensamer R. Compressed sensing in astronomy,. IEEE Journal of Selected Topics in Signal Processing 2009;28(2):609–18. [8] Duarte M, Davenport M, Takhar D, Laska J, Sun T, Kelly K, et al. Single-pixel imaging via compressive sampling,. IEEE Signal Processing Magazine 2008;25(2):83–91. [9] Han B, Zhou CL. High-resolution infrared or millimeter wave imaging based on compressed sensing. Infrared 2010;31(8):1–5 in Chinese. [10] Marcia Roummel F, Harmany Zachary T, Willett Rebecca M. Compressive coded aperture imaging. In: Proceeding of SPIE; 2009. p. 72460G. [11] Willett Rebecca M, Marcia Roummel F, Nichols Jonathan M. Compressed sensing for practical optical imaging systems: a tutorial. Optical Engineering 2011;50(7):072601. [12] Candes E. The restricted isometry property and its implications for compressed sensing. Compte Rendus de l’Academie des Sciences 2006;346(I): 589–92. [13] Candes E, Tao T. Decoding by linear programming,. IEEE Transaction on Information Theory 2005;51(12):4203–15. [14] Candes E, Romberg J, Tao T. Stable signal recovery from incomplete and inaccurate measurements. Communications on Pure and Applied Mathematics 2006;59(8):1207–23. [15] Donoho DL. For most large underdetermined systems of linear equations, the minimal L1 norm solution is also the sparsest solution. Communications On Pure And Applied Mathematics 2006;59(6):797–829. [16] Tropp J, Gilbert A. Signal recovery from random measurements via orthogonal matching pursuit. Transactions on Information Theory 2007;53(12):4655–66. [17] Herrity, K, Gilbert, A and Tropp, JA. Sparse approximation via iterative thresholding, in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, Washington DC, USA, 2006, 624-627. [18] Candes E, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,. IEEE Transactions on Information Theory 2006;52(2):489–509. [19] Portnoy AD, Pitsianis NP, Brady DJ, Guo J, Fiddy MA, Feldman MR, et al. Thin digital imaging system using focal plane coding. Proc. of SPIE-IS&T Electronic Imaging, SPIE 2006;6065:108–15. [20] Figueiredo MA, Nowak R, Wright SJ. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems, IEEE. Journal of Selected Topics in Signal Processing 2007;1(4):586–98. [21] Voelz David. Computational Fourier Optics_A Matlab Tutorial. Washington: SPIE Press; 2010.