6466
OPTICS LETTERS / Vol. 39, No. 22 / November 15, 2014
Single-shot phase imaging with a coded aperture Ryoichi Horisaki,* Yusuke Ogura, Masahiko Aino, and Jun Tanida Department of Information and Physical Sciences, Graduate School of Information Science and Technology, Osaka University, 1-5 Yamadaoka, Suita, Osaka 565-0871, Japan *Corresponding author: [email protected]‑u.ac.jp Received August 7, 2014; revised September 26, 2014; accepted October 12, 2014; posted October 14, 2014 (Doc. ID 220578); published November 11, 2014 We present a method of quantitatively acquiring a large complex field, containing not only amplitude information but also phase information, based on single-shot phase imaging with a coded aperture (SPICA). In SPICA, the propagating field from an object illuminated by partially coherent visible light is sieved by a coded mask, and the sieved field propagates to an image sensor, where it is captured. The sieved field is recovered from the single captured intensity image via a phase retrieval algorithm with an amplitude support constraint using the mask pattern, and then the object’s complex field is reconstructed from the recovered sieved field by an algorithm employing a sparsity constraint based on compressive sensing. The system model and the theoretical bounds of SPICA are derived. We also verified the concept with numerical demonstrations. © 2014 Optical Society of America OCIS codes: (110.1758) Computational imaging; (170.1630) Coded aperture imaging; (100.5070) Phase retrieval; (090.1995) Digital holography. http://dx.doi.org/10.1364/OL.39.006466
A light wave is represented by amplitude and phase, and together these are treated as a complex field [1]. The complex field from an object is generally observed using interferometric approaches, such as holography. Many interferometric approaches have been proposed and have realized various useful functions, including threedimensional imaging and phase imaging [2,3]. They are especially useful for biomedical imaging because most of the objects in this field are three-dimensional and transparent. In general, interferometric systems require a reference optical path or multiple shots with a phaseshift to determine the interference quantitatively, but this increases the size and complexity of the hardware or the observation time. An exception is Gabor holography, because in this technique, the light that illuminates the object and the reference light are shared [4]. However, it suffers from autocorrelation and conjugation of the object field. Another exception is coherent diffraction imaging (CDI), where the complex field is recovered from an intensity image of a diffraction pattern of an object illuminated by coherent light by using an algorithm for phase retrieval (PR) [5,6]. However, when recovering the complex field, in order to constrain the field in the PR process, the object must be much smaller than the image sensor used, or multishot imaging must be executed with different imaging conditions, which is achieved, for example, by scanning the field-of-view or the focusing distance, by sequential modulation of the amplitude or phase, or by using multiple structured illuminations [7–12]. CDI has been mainly used in the X-ray region, where lenses are difficult to fabricate. Here, we propose single-shot CDI with a coded aperture (CA). We refer to the method as single-shot phase imaging with a coded aperture (SPICA). The CA consists of many specially or randomly arranged pinholes in a mask plate [13]. CA imaging has been typically used instead of a pinhole camera to increase the light efficiency, especially in the X-ray region. In the visible region, however, CA imaging suffers from diffraction by the pinholes. Thus, multishot imaging with different CAs has been proposed and demonstrated in the visible region [14]. 0146-9592/14/226466-04$15.00/0
Regularly arranged pinholes have been applied to CDI to enable use of a noniterative PR process for a small object [7,10]. SPICA realizes single-shot imaging for a complex field by using only a CA, an image sensor, and a partially coherent light source in the visible region. Furthermore, as demonstrated in this Letter, the acquirable complex field in SPICA is the same size as the image sensor and is larger than that of conventional single-shot CDI. Thus, SPICA significantly extends the imaging capability of CDI with a low-cost approach. SPICA consists of compressive Fresnel holography (CFH) and CDI with a CA, as shown in Fig. 1. CFH is a holographic imaging technique in which the propagating field from an object is sparsely sampled with randomly arranged detectors to reduce redundancy in the measurements, and the object is reconstructed from the sparse sampling data with a numerical algorithm employing the sparsity of the object based on compressive sensing (CS) [15]. CS is a framework for observing a large-sized object with fewer measurements compared with the conventional sampling theorem [16]. In SPICA, the complex field F of the object propagates to the CA plane. The propagating field is multiplied by the mask Compressive Fresnel holography
Coherent diffraction imaging with a coded apeture (CA)
Reconstruction with Phase retrieval with sparsity of the object CA-based support
Fig. 1. Schematic diagram of SPICA. The object field (s–t plane) propagates to the CA (u–v plane) and is sieved by the CA. The sieved field propagates to the sensor (x–y plane). © 2014 Optical Society of America
November 15, 2014 / Vol. 39, No. 22 / OPTICS LETTERS
pattern M of the CA, thus sieving the field. The sieved complex field A just behind the CA propagates to the image sensor, which captures a single intensity image G. The former part corresponds to CFH, and the latter corresponds to CDI. As shown in Fig. 1, the sieved field A is recovered from the single intensity image G by using a PR algorithm. The object field F is reconstructed from the recovered sieved field with a CS algorithm. In CS, the object must be sparse in a certain domain, e.g., real space and wavelet space, but speckles on the object reduce the sparsity. To reduce the speckles and increase the sparsity, we assume partially coherent light, e.g., from a light-emitting diode (LED), to illuminate the object. The t, v, and y-axes of the object, the CA, and the sensor planes in Fig. 1 are omitted for simplicity. Then the system model of SPICA is written as Z Au Mu
Pu − s; zF ; λc F sds;
2 Z Gx Px − u; zA ; λc Audu ;
(1)
(2)
where P•; z; λ is the point-spread function (PSF) of the Fresnel propagation at distance z and wavelength λ [1], zF is the distance between the object and the CA, zA is the distance between the CA and the image sensor, and s, u, and x are the spatial axes of the object, the CA, and the sensor planes, as shown in Fig. 1. We assume an approximately spatially coherent quasi-monochromatic light source with a center wavelength λc to illuminate the object. This system model can be readily extended to multidimensional imaging, including three spatial dimensions and the wavelength dimension [17]. Equation (1) is rewritten in matrix notation as a MPzF ;λc f ;
(3)
6467
(1) p An initial complex field is defined as the amplitude g with a random phase distribution. (2) The complex field on the sensor plane is back-propagated by distance zA . (3) The back-propagated field is multiplied with M, and the result is defined as the new estimation aˆ . (4) aˆ is forward-propagated onto the sensor plane, pand the amplitude of the resultant field is replaced by g. Steps (2) to (4) are iterated until aˆ converges. In the PR, the support constraint is flexible, and it also can constrain the complex field on the sensor plane. Such a support constraint is useful for treating the extended propagating field on the sensor plane. The inversion of Eq. (3) is ill-posed because the mask pattern M contains zeros. To solve the inversion problem, we use a CS algorithm called two-step iterative shrinkage/thresholding (TwIST) [18]. TwIST is an iterative convex optimization algorithm using two previous estimates to improve the convergence in the iterations. TwIST reconstructs the original field f as fˆ by solving the following optimization problem: fˆ argmin∥ˆa − MPzF ;λc f ∥l2 τRf ; f
(5)
where ∥ • ∥l2 is the l2 norm, τ is a regularization parameter, and R• is a regularizer. SPICA was numerically implemented, and the results are shown in Fig. 2. The pixel count N 2x was 5122 . The pixel pitch δx of the object, CA, and sensor planes was 5 μm. The wavelength λc was 532 nm. The amplitude and phase of the object field f are shown in Figs. 2(a) and 2(b), respectively. These were composed of a standard natural image: Pepper. The CA with the mask M in Fig. 2(c) was located 50 mm (zF ) from the object. The number N p of white pixels in Fig. 2(c), which are the randomly arranged pinholes in the mask where light passes through, was 65536, which is 25% of all pixels (5122 ). The image sensor was located 1 mm (zA ) from the CA.
where a ∈ CN x ×1 is the vectored sieved complex field just behind the CA, f ∈ CN x ×1 is the vectored object field, N x is the number of elements of the one-dimensional data, M ∈ RN x ×N x is a diagonal matrix in which the diagonal elements correspond to the mask pattern, and Pz;λ ∈ CN x ×N x is a Toeplitz matrix expressing the convolution of the PSF of the Fresnel propagation at distance z and wavelength λ. Equation (2) can be rewritten as g jPzA ;λc aj2 ;
(4)
where g ∈ RN x ×1 is the vectored captured intensity data. The inversion of Eq. (4) can be solved by using the input-output PR (IOPR) algorithm [6]. IOPR iteratively retrieves the object’s complex field from its propagated intensity image with a support constraint of the object’s amplitude. In practice, however, the object must be very small to satisfy this constraint. In our case, the sieved field a is spatially sparse, and available points on the field are known because of the known mask pattern M of the CA. Thus, our method alleviates the limitation regarding the object’s size in conventional CDI. IOPR with the CA-based support constraint is executed as follows:
Fig. 2. Simulation results. The (a) amplitude and (b) phase of the object field. (c) The mask pattern of the CA, where the upper-right part shows a magnified 10 × 10 pattern on the CA. (d) The captured intensity image. The (e) amplitude and (f) phase of the reconstructed object field. Phases are normalized in the interval −π; π.
6468
OPTICS LETTERS / Vol. 39, No. 22 / November 15, 2014
The captured intensity image g is shown in Fig. 2(d), where Gaussian random noise with a signal-to-noise ratio (SNR) of 30 dB was added. The amplitude and phase after phase bias matching of the object field fˆ reconstructed by IOPR and TwIST are shown in Figs. 2(e) and 2(f), respectively. We chose the two-dimensional total variation (2DTV) as the regularizer R in TwIST [19]. By comparing Figs. 2(a) and 2(e), and Figs. 2(b) and 2(f), the amplitude and the relative phase were recovered successfully. The peak SNR (PSNR) between the original and reconstructed object fields was 29.7 dB. The relationship between the measurement noise level (SNR) and the reconstruction fidelity (PNSR) in the simulated setup is plotted in Fig. 3. The plot shows averaged PSNRs obtained with ten different mask patterns. The PSNRs rise at the measurement noise level of around 20 dB. The theoretical bounds of SPICA were determined by combining those of the inversion of CFH in Eq. (3) and the inversion of CDI in Eq. (4). The following discussions of the bounds assume an ideal noiseless measurement. The upper bound of CFH has been derived in [20]. When the distance zF between the object and the CA plane satisfies the condition zF ≥
N x δ2x ; λc
(6)
the number N CFH of required sampling points in CFH is given by N CFH ≥ ck log N 2x ;
(7)
Reconstruction fidelity: PSNR (dB)
where c is a constant depending on the acceptable reconstruction fidelity, and k shows the sparsity, which is the number of nonzero coefficients in the regularizer domain. In the case of SPICA, N p N CFH . The maximum numerical aperture (NA) is introduced with the minimum distance zF in Eq. (6), and it is given by N x δx ∕2zF λc ∕2δx . Then the achievable resolution is calculated from 0.5λc ∕NA δx . A point source with a size of δx on the object plane is resolved in this case, and this resolution is equivalent to that in conventional imaging systems [1]. Additional optics also can be employed to magnify the object optically. In the first simulation in 35 30 25 20 15 10 5 0
0
20 40 60 Measurement noise level: SNR (dB)
Fig. 3. Relationship between the measurement noise level and the reconstruction fidelity.
Fig. 2, the right side of Eq. (6) became 24.1 mm, and the condition in Eq. (6) was satisfied. The upper bound of PR in CDI was shown in [21]. To recover complex data from the magnitude of its Fourier transform, oversampling is necessary to find the unique autocorrelation of the object. When the oversampling factor is defined as ρ, the condition for the distance zA between the CA and sensor planes is zA ≥
p 2 ρδx : λc
(8)
In this case, a single pixel on the CA contributes more than ρ pixels on the image sensor. The condition for the number N p of pinholes in the CA in SPICA is given by Np ≤
N 2x : ρ
(9)
In [21], twice-as-large oversampling in each dimension is required because, in each dimension, the size of the autocorrelation is twice as large as that of the object. Since, the CA is two-dimensional, ρ ≈ 4. The discussions in this paragraph enable us to decide the distance zA without relying on the object and the CFH process. It indicates the possibility of reducing the thickness of the hardware, approximating the size of the propagating field on the sensor to the size of the CA, and increasing the SNR in the measurement. The right sides of Eqs. (8) and (9) were 94.0 μm and 65536 in the case of the first simulation, and thus the conditions were satisfied. By combining Eqs. (7) and (9), the upper bound of SPICA is written as 1 ck ≥ cq; ρ log N 2x N 2x
(10)
where q k∕N 2x is defined as the relative sparsity. The relationship between the relative sparsity and the reconstruction fidelity (PNSR) for different oversampling factors ρ was simulated, and the results are plotted in Fig. 4. Here, the pixel count N 2x was 1282 . The number N p of pinholes was set as the upper bound of Eq. (9). The object in this simulation consisted of randomly distributed point sources having the same intensity and different random phases. The other parameters were the same as those in the first simulation. The distances zF and zA satisfied the conditions in Eqs. (6) and (8) in this case. We chose the l1 norm in real space as the regularizer R in TwIST because the object was spatially sparse. Thus, the sparsity k is equal to the number of point sources. Ten distributions of point sources and ten mask patterns were realized at each relative sparsity q, and the averages were plotted. The reconstruction fidelity was evaluated using the PSNR between the original and reconstructed fields. Figure 4 shows that the PSNR for each oversampling factor ρ falls sharply at different relative sparsities q. If we assume a PSNR of 40 dB as the acceptable reconstruction fidelity, the critical relative sparsity q and the constant c are as shown in Table 1, calculated from Fig. 4 and Eq. (10). The constant c has similar values at the different oversampling factors ρ. This supports the
Reconstruction fidelity: PSNR (dB)
November 15, 2014 / Vol. 39, No. 22 / OPTICS LETTERS 70
ρ=4 ρ=5 ρ=6 ρ=7 ρ=8
60 50 40 30 20 10 0 0
0.05 0.1 Relative sparsity
0.15
Fig. 4. Relationship between the relative sparsity and the reconstruction fidelity.
Table 1. Calculated Bounds of SPICA ρ
4
5
6
7
8
q c
0.063 0.41
0.048 0.43
0.038 0.45
0.032 0.46
0.028 0.46
derivation of the theoretical bounds. One approach to increase the critical relative sparsity q is to reduce the oversampling factor ρ. In this Letter, although we used the IOPR, which is a very basic PR algorithm, it should be possible to use state-of-the-art PR methods for SPICA and to reduce the oversampling factor ρ. In this Letter, we proposed single-shot coherent diffraction imaging (CDI) with a coded aperture (CA), which is referred to as single-shot phase imaging with a coded aperture (SPICA). SPICA combines compressive Fresnel holography (CFH) and CDI by using the CA. The optical setup of SPICA is composed of only the CA, an image sensor, and a partially coherent light source. An object is illuminated by the light source to reduce the speckles and increase the sparsity of the object in a certain domain. The object field propagates to the image sensor via the CFH process with the CA, and the intensity of a single diffraction pattern is captured. The sieved complex field just behind the CA is recovered from the single intensity image by a phase retrieval (PR) algorithm with a CA-based support constraint. The original field is reconstructed from the recovered sieved complex field by a sparsity-based algorithm used in compressive sensing (CS). We derived the system model and the theoretical bounds of SPICA. The concept of SPICA was verified by simulations. A future issue is to consider more-sophisticated analysis of the performance of SPICA. For example,
6469
the robustness against noise and the stability of convergence in the reconstruction are important in practical situations. State-of-the-art algorithms for CS and PR will be adapted to SPICA to alleviate the theoretical bounds discussed here. We also will investigate optimization of the CA and imaging conditions. SPICA realized single-shot quantitative imaging for a large complex field with a compact and simple setup. Although two-dimensional imaging was demonstrated in this Letter, SPICA is extendable to higher dimensional imaging, including the depth and wavelength dimensions. SPICA is also suitable for fluorescence imaging and for CDI in various spectral regions, such as X-ray and electron imaging. Thus, SPICA is a promising new platform for future imaging systems. References 1. J. W. Goodman, Introduction to Fourier Optics (McGrawHill, 1996). 2. G. Nehmetallah and P. P. Banerjee, Adv. Opt. Photon. 4, 472 (2012). 3. B. Bhaduri, C. Edwards, H. Pham, R. Zhou, T. H. Nguyen, L. L. Goddard, and G. Popescu, Adv. Opt. Photon. 6, 57 (2014). 4. D. Gabor, Nature 161, 777 (1948). 5. S. Marchesini, H. Chapman, S. Hau-Riege, R. London, A. Szoke, H. He, M. Howells, H. Padmore, R. Rosen, J. Spence, and U. Weierstall, Opt. Express 11, 2344 (2003). 6. J. R. Fienup, Appl. Opt. 21, 2758 (1982). 7. N. Nakajima, Phys. Rev. Lett. 98, 223901 (2007). 8. P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, Science 321, 379 (2008). 9. N. Streibl, Opt. Commun. 49, 6 (1984). 10. G.-X. Wei, L.-L. Lu, C.-S. Guo, and H.-T. Wang, Appl. Opt. 48, 5099 (2009). 11. C. Falldorf, M. Agour, C. v. Kopylow, and R. B. Bergmann, Appl. Opt. 49, 1826 (2010). 12. A. Fannjiang and W. Liao, J. Opt. Soc. Am. A 29, 1847 (2012). 13. E. M. Fenimore and T. M. Cannon, Appl. Opt. 17, 337 (1978). 14. A. Busboom, H. D. Schotten, and H. Elders-Boll, J. Opt. Soc. Am. A 14, 1058 (1997). 15. Y. Rivenson, A. Stern, and B. Javidi, J. Disp. Technol. 6, 506 (2010). 16. D. L. Donoho, IEEE Trans. Inf. Theory 52, 1289 (2006). 17. R. Horisaki, J. Tanida, A. Stern, and B. Javidi, Opt. Lett. 37, 2013 (2012). 18. J. M. Bioucas-Dias and M. A. T. Figueiredo, IEEE Trans. Image Process. 16, 2992 (2007). 19. L. I. Rudin, S. Osher, and E. Fatemi, Physica D 60, 259 (1992). 20. Y. Rivenson and A. Stern, Opt. Lett. 36, 3365 (2011). 21. R. H. T. Bates, Comput. Vis. Graph. Image Process. 25, 205 (1984).
OPTICS LETTERS / Vol. 39, No. 22 / November 15, 2014
Single-shot phase imaging with a coded aperture Ryoichi Horisaki,* Yusuke Ogura, Masahiko Aino, and Jun Tanida Department of Information and Physical Sciences, Graduate School of Information Science and Technology, Osaka University, 1-5 Yamadaoka, Suita, Osaka 565-0871, Japan *Corresponding author: [email protected]‑u.ac.jp Received August 7, 2014; revised September 26, 2014; accepted October 12, 2014; posted October 14, 2014 (Doc. ID 220578); published November 11, 2014 We present a method of quantitatively acquiring a large complex field, containing not only amplitude information but also phase information, based on single-shot phase imaging with a coded aperture (SPICA). In SPICA, the propagating field from an object illuminated by partially coherent visible light is sieved by a coded mask, and the sieved field propagates to an image sensor, where it is captured. The sieved field is recovered from the single captured intensity image via a phase retrieval algorithm with an amplitude support constraint using the mask pattern, and then the object’s complex field is reconstructed from the recovered sieved field by an algorithm employing a sparsity constraint based on compressive sensing. The system model and the theoretical bounds of SPICA are derived. We also verified the concept with numerical demonstrations. © 2014 Optical Society of America OCIS codes: (110.1758) Computational imaging; (170.1630) Coded aperture imaging; (100.5070) Phase retrieval; (090.1995) Digital holography. http://dx.doi.org/10.1364/OL.39.006466
A light wave is represented by amplitude and phase, and together these are treated as a complex field [1]. The complex field from an object is generally observed using interferometric approaches, such as holography. Many interferometric approaches have been proposed and have realized various useful functions, including threedimensional imaging and phase imaging [2,3]. They are especially useful for biomedical imaging because most of the objects in this field are three-dimensional and transparent. In general, interferometric systems require a reference optical path or multiple shots with a phaseshift to determine the interference quantitatively, but this increases the size and complexity of the hardware or the observation time. An exception is Gabor holography, because in this technique, the light that illuminates the object and the reference light are shared [4]. However, it suffers from autocorrelation and conjugation of the object field. Another exception is coherent diffraction imaging (CDI), where the complex field is recovered from an intensity image of a diffraction pattern of an object illuminated by coherent light by using an algorithm for phase retrieval (PR) [5,6]. However, when recovering the complex field, in order to constrain the field in the PR process, the object must be much smaller than the image sensor used, or multishot imaging must be executed with different imaging conditions, which is achieved, for example, by scanning the field-of-view or the focusing distance, by sequential modulation of the amplitude or phase, or by using multiple structured illuminations [7–12]. CDI has been mainly used in the X-ray region, where lenses are difficult to fabricate. Here, we propose single-shot CDI with a coded aperture (CA). We refer to the method as single-shot phase imaging with a coded aperture (SPICA). The CA consists of many specially or randomly arranged pinholes in a mask plate [13]. CA imaging has been typically used instead of a pinhole camera to increase the light efficiency, especially in the X-ray region. In the visible region, however, CA imaging suffers from diffraction by the pinholes. Thus, multishot imaging with different CAs has been proposed and demonstrated in the visible region [14]. 0146-9592/14/226466-04$15.00/0
Regularly arranged pinholes have been applied to CDI to enable use of a noniterative PR process for a small object [7,10]. SPICA realizes single-shot imaging for a complex field by using only a CA, an image sensor, and a partially coherent light source in the visible region. Furthermore, as demonstrated in this Letter, the acquirable complex field in SPICA is the same size as the image sensor and is larger than that of conventional single-shot CDI. Thus, SPICA significantly extends the imaging capability of CDI with a low-cost approach. SPICA consists of compressive Fresnel holography (CFH) and CDI with a CA, as shown in Fig. 1. CFH is a holographic imaging technique in which the propagating field from an object is sparsely sampled with randomly arranged detectors to reduce redundancy in the measurements, and the object is reconstructed from the sparse sampling data with a numerical algorithm employing the sparsity of the object based on compressive sensing (CS) [15]. CS is a framework for observing a large-sized object with fewer measurements compared with the conventional sampling theorem [16]. In SPICA, the complex field F of the object propagates to the CA plane. The propagating field is multiplied by the mask Compressive Fresnel holography
Coherent diffraction imaging with a coded apeture (CA)
Reconstruction with Phase retrieval with sparsity of the object CA-based support
Fig. 1. Schematic diagram of SPICA. The object field (s–t plane) propagates to the CA (u–v plane) and is sieved by the CA. The sieved field propagates to the sensor (x–y plane). © 2014 Optical Society of America
November 15, 2014 / Vol. 39, No. 22 / OPTICS LETTERS
pattern M of the CA, thus sieving the field. The sieved complex field A just behind the CA propagates to the image sensor, which captures a single intensity image G. The former part corresponds to CFH, and the latter corresponds to CDI. As shown in Fig. 1, the sieved field A is recovered from the single intensity image G by using a PR algorithm. The object field F is reconstructed from the recovered sieved field with a CS algorithm. In CS, the object must be sparse in a certain domain, e.g., real space and wavelet space, but speckles on the object reduce the sparsity. To reduce the speckles and increase the sparsity, we assume partially coherent light, e.g., from a light-emitting diode (LED), to illuminate the object. The t, v, and y-axes of the object, the CA, and the sensor planes in Fig. 1 are omitted for simplicity. Then the system model of SPICA is written as Z Au Mu
Pu − s; zF ; λc F sds;
2 Z Gx Px − u; zA ; λc Audu ;
(1)
(2)
where P•; z; λ is the point-spread function (PSF) of the Fresnel propagation at distance z and wavelength λ [1], zF is the distance between the object and the CA, zA is the distance between the CA and the image sensor, and s, u, and x are the spatial axes of the object, the CA, and the sensor planes, as shown in Fig. 1. We assume an approximately spatially coherent quasi-monochromatic light source with a center wavelength λc to illuminate the object. This system model can be readily extended to multidimensional imaging, including three spatial dimensions and the wavelength dimension [17]. Equation (1) is rewritten in matrix notation as a MPzF ;λc f ;
(3)
6467
(1) p An initial complex field is defined as the amplitude g with a random phase distribution. (2) The complex field on the sensor plane is back-propagated by distance zA . (3) The back-propagated field is multiplied with M, and the result is defined as the new estimation aˆ . (4) aˆ is forward-propagated onto the sensor plane, pand the amplitude of the resultant field is replaced by g. Steps (2) to (4) are iterated until aˆ converges. In the PR, the support constraint is flexible, and it also can constrain the complex field on the sensor plane. Such a support constraint is useful for treating the extended propagating field on the sensor plane. The inversion of Eq. (3) is ill-posed because the mask pattern M contains zeros. To solve the inversion problem, we use a CS algorithm called two-step iterative shrinkage/thresholding (TwIST) [18]. TwIST is an iterative convex optimization algorithm using two previous estimates to improve the convergence in the iterations. TwIST reconstructs the original field f as fˆ by solving the following optimization problem: fˆ argmin∥ˆa − MPzF ;λc f ∥l2 τRf ; f
(5)
where ∥ • ∥l2 is the l2 norm, τ is a regularization parameter, and R• is a regularizer. SPICA was numerically implemented, and the results are shown in Fig. 2. The pixel count N 2x was 5122 . The pixel pitch δx of the object, CA, and sensor planes was 5 μm. The wavelength λc was 532 nm. The amplitude and phase of the object field f are shown in Figs. 2(a) and 2(b), respectively. These were composed of a standard natural image: Pepper. The CA with the mask M in Fig. 2(c) was located 50 mm (zF ) from the object. The number N p of white pixels in Fig. 2(c), which are the randomly arranged pinholes in the mask where light passes through, was 65536, which is 25% of all pixels (5122 ). The image sensor was located 1 mm (zA ) from the CA.
where a ∈ CN x ×1 is the vectored sieved complex field just behind the CA, f ∈ CN x ×1 is the vectored object field, N x is the number of elements of the one-dimensional data, M ∈ RN x ×N x is a diagonal matrix in which the diagonal elements correspond to the mask pattern, and Pz;λ ∈ CN x ×N x is a Toeplitz matrix expressing the convolution of the PSF of the Fresnel propagation at distance z and wavelength λ. Equation (2) can be rewritten as g jPzA ;λc aj2 ;
(4)
where g ∈ RN x ×1 is the vectored captured intensity data. The inversion of Eq. (4) can be solved by using the input-output PR (IOPR) algorithm [6]. IOPR iteratively retrieves the object’s complex field from its propagated intensity image with a support constraint of the object’s amplitude. In practice, however, the object must be very small to satisfy this constraint. In our case, the sieved field a is spatially sparse, and available points on the field are known because of the known mask pattern M of the CA. Thus, our method alleviates the limitation regarding the object’s size in conventional CDI. IOPR with the CA-based support constraint is executed as follows:
Fig. 2. Simulation results. The (a) amplitude and (b) phase of the object field. (c) The mask pattern of the CA, where the upper-right part shows a magnified 10 × 10 pattern on the CA. (d) The captured intensity image. The (e) amplitude and (f) phase of the reconstructed object field. Phases are normalized in the interval −π; π.
6468
OPTICS LETTERS / Vol. 39, No. 22 / November 15, 2014
The captured intensity image g is shown in Fig. 2(d), where Gaussian random noise with a signal-to-noise ratio (SNR) of 30 dB was added. The amplitude and phase after phase bias matching of the object field fˆ reconstructed by IOPR and TwIST are shown in Figs. 2(e) and 2(f), respectively. We chose the two-dimensional total variation (2DTV) as the regularizer R in TwIST [19]. By comparing Figs. 2(a) and 2(e), and Figs. 2(b) and 2(f), the amplitude and the relative phase were recovered successfully. The peak SNR (PSNR) between the original and reconstructed object fields was 29.7 dB. The relationship between the measurement noise level (SNR) and the reconstruction fidelity (PNSR) in the simulated setup is plotted in Fig. 3. The plot shows averaged PSNRs obtained with ten different mask patterns. The PSNRs rise at the measurement noise level of around 20 dB. The theoretical bounds of SPICA were determined by combining those of the inversion of CFH in Eq. (3) and the inversion of CDI in Eq. (4). The following discussions of the bounds assume an ideal noiseless measurement. The upper bound of CFH has been derived in [20]. When the distance zF between the object and the CA plane satisfies the condition zF ≥
N x δ2x ; λc
(6)
the number N CFH of required sampling points in CFH is given by N CFH ≥ ck log N 2x ;
(7)
Reconstruction fidelity: PSNR (dB)
where c is a constant depending on the acceptable reconstruction fidelity, and k shows the sparsity, which is the number of nonzero coefficients in the regularizer domain. In the case of SPICA, N p N CFH . The maximum numerical aperture (NA) is introduced with the minimum distance zF in Eq. (6), and it is given by N x δx ∕2zF λc ∕2δx . Then the achievable resolution is calculated from 0.5λc ∕NA δx . A point source with a size of δx on the object plane is resolved in this case, and this resolution is equivalent to that in conventional imaging systems [1]. Additional optics also can be employed to magnify the object optically. In the first simulation in 35 30 25 20 15 10 5 0
0
20 40 60 Measurement noise level: SNR (dB)
Fig. 3. Relationship between the measurement noise level and the reconstruction fidelity.
Fig. 2, the right side of Eq. (6) became 24.1 mm, and the condition in Eq. (6) was satisfied. The upper bound of PR in CDI was shown in [21]. To recover complex data from the magnitude of its Fourier transform, oversampling is necessary to find the unique autocorrelation of the object. When the oversampling factor is defined as ρ, the condition for the distance zA between the CA and sensor planes is zA ≥
p 2 ρδx : λc
(8)
In this case, a single pixel on the CA contributes more than ρ pixels on the image sensor. The condition for the number N p of pinholes in the CA in SPICA is given by Np ≤
N 2x : ρ
(9)
In [21], twice-as-large oversampling in each dimension is required because, in each dimension, the size of the autocorrelation is twice as large as that of the object. Since, the CA is two-dimensional, ρ ≈ 4. The discussions in this paragraph enable us to decide the distance zA without relying on the object and the CFH process. It indicates the possibility of reducing the thickness of the hardware, approximating the size of the propagating field on the sensor to the size of the CA, and increasing the SNR in the measurement. The right sides of Eqs. (8) and (9) were 94.0 μm and 65536 in the case of the first simulation, and thus the conditions were satisfied. By combining Eqs. (7) and (9), the upper bound of SPICA is written as 1 ck ≥ cq; ρ log N 2x N 2x
(10)
where q k∕N 2x is defined as the relative sparsity. The relationship between the relative sparsity and the reconstruction fidelity (PNSR) for different oversampling factors ρ was simulated, and the results are plotted in Fig. 4. Here, the pixel count N 2x was 1282 . The number N p of pinholes was set as the upper bound of Eq. (9). The object in this simulation consisted of randomly distributed point sources having the same intensity and different random phases. The other parameters were the same as those in the first simulation. The distances zF and zA satisfied the conditions in Eqs. (6) and (8) in this case. We chose the l1 norm in real space as the regularizer R in TwIST because the object was spatially sparse. Thus, the sparsity k is equal to the number of point sources. Ten distributions of point sources and ten mask patterns were realized at each relative sparsity q, and the averages were plotted. The reconstruction fidelity was evaluated using the PSNR between the original and reconstructed fields. Figure 4 shows that the PSNR for each oversampling factor ρ falls sharply at different relative sparsities q. If we assume a PSNR of 40 dB as the acceptable reconstruction fidelity, the critical relative sparsity q and the constant c are as shown in Table 1, calculated from Fig. 4 and Eq. (10). The constant c has similar values at the different oversampling factors ρ. This supports the
Reconstruction fidelity: PSNR (dB)
November 15, 2014 / Vol. 39, No. 22 / OPTICS LETTERS 70
ρ=4 ρ=5 ρ=6 ρ=7 ρ=8
60 50 40 30 20 10 0 0
0.05 0.1 Relative sparsity
0.15
Fig. 4. Relationship between the relative sparsity and the reconstruction fidelity.
Table 1. Calculated Bounds of SPICA ρ
4
5
6
7
8
q c
0.063 0.41
0.048 0.43
0.038 0.45
0.032 0.46
0.028 0.46
derivation of the theoretical bounds. One approach to increase the critical relative sparsity q is to reduce the oversampling factor ρ. In this Letter, although we used the IOPR, which is a very basic PR algorithm, it should be possible to use state-of-the-art PR methods for SPICA and to reduce the oversampling factor ρ. In this Letter, we proposed single-shot coherent diffraction imaging (CDI) with a coded aperture (CA), which is referred to as single-shot phase imaging with a coded aperture (SPICA). SPICA combines compressive Fresnel holography (CFH) and CDI by using the CA. The optical setup of SPICA is composed of only the CA, an image sensor, and a partially coherent light source. An object is illuminated by the light source to reduce the speckles and increase the sparsity of the object in a certain domain. The object field propagates to the image sensor via the CFH process with the CA, and the intensity of a single diffraction pattern is captured. The sieved complex field just behind the CA is recovered from the single intensity image by a phase retrieval (PR) algorithm with a CA-based support constraint. The original field is reconstructed from the recovered sieved complex field by a sparsity-based algorithm used in compressive sensing (CS). We derived the system model and the theoretical bounds of SPICA. The concept of SPICA was verified by simulations. A future issue is to consider more-sophisticated analysis of the performance of SPICA. For example,
6469
the robustness against noise and the stability of convergence in the reconstruction are important in practical situations. State-of-the-art algorithms for CS and PR will be adapted to SPICA to alleviate the theoretical bounds discussed here. We also will investigate optimization of the CA and imaging conditions. SPICA realized single-shot quantitative imaging for a large complex field with a compact and simple setup. Although two-dimensional imaging was demonstrated in this Letter, SPICA is extendable to higher dimensional imaging, including the depth and wavelength dimensions. SPICA is also suitable for fluorescence imaging and for CDI in various spectral regions, such as X-ray and electron imaging. Thus, SPICA is a promising new platform for future imaging systems. References 1. J. W. Goodman, Introduction to Fourier Optics (McGrawHill, 1996). 2. G. Nehmetallah and P. P. Banerjee, Adv. Opt. Photon. 4, 472 (2012). 3. B. Bhaduri, C. Edwards, H. Pham, R. Zhou, T. H. Nguyen, L. L. Goddard, and G. Popescu, Adv. Opt. Photon. 6, 57 (2014). 4. D. Gabor, Nature 161, 777 (1948). 5. S. Marchesini, H. Chapman, S. Hau-Riege, R. London, A. Szoke, H. He, M. Howells, H. Padmore, R. Rosen, J. Spence, and U. Weierstall, Opt. Express 11, 2344 (2003). 6. J. R. Fienup, Appl. Opt. 21, 2758 (1982). 7. N. Nakajima, Phys. Rev. Lett. 98, 223901 (2007). 8. P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, Science 321, 379 (2008). 9. N. Streibl, Opt. Commun. 49, 6 (1984). 10. G.-X. Wei, L.-L. Lu, C.-S. Guo, and H.-T. Wang, Appl. Opt. 48, 5099 (2009). 11. C. Falldorf, M. Agour, C. v. Kopylow, and R. B. Bergmann, Appl. Opt. 49, 1826 (2010). 12. A. Fannjiang and W. Liao, J. Opt. Soc. Am. A 29, 1847 (2012). 13. E. M. Fenimore and T. M. Cannon, Appl. Opt. 17, 337 (1978). 14. A. Busboom, H. D. Schotten, and H. Elders-Boll, J. Opt. Soc. Am. A 14, 1058 (1997). 15. Y. Rivenson, A. Stern, and B. Javidi, J. Disp. Technol. 6, 506 (2010). 16. D. L. Donoho, IEEE Trans. Inf. Theory 52, 1289 (2006). 17. R. Horisaki, J. Tanida, A. Stern, and B. Javidi, Opt. Lett. 37, 2013 (2012). 18. J. M. Bioucas-Dias and M. A. T. Figueiredo, IEEE Trans. Image Process. 16, 2992 (2007). 19. L. I. Rudin, S. Osher, and E. Fatemi, Physica D 60, 259 (1992). 20. Y. Rivenson and A. Stern, Opt. Lett. 36, 3365 (2011). 21. R. H. T. Bates, Comput. Vis. Graph. Image Process. 25, 205 (1984).
