Multidimensional object acquisition by single-shot phase imaging with a coded aperture Ryoichi Horisaki∗ 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 ∗ [email protected]
Abstract: We propose a generalized framework for quantitatively acquiring multidimensional complex objects based on single-shot phase imaging with a coded aperture (SPICA). In multidimensional SPICA, a propagating field from a multidimensional complex object is sieved by a coded aperture, the sieved field is modulated by an optical element, which is called coding optics, and then the resultant field is captured by a monochrome image sensor. The original complex field is reconstructed from the single captured intensity image by a phase retrieval algorithm with a support constraint of the coded aperture and a sparsity-based reconstruction algorithm based on compressive sensing. We also present theoretical conditions for the proposed method. As a demonstration, we numerically verified an application of this generalized framework for single-shot acquisition of depth-variant multispectral objects. © 2015 Optical Society of America OCIS codes: (110.1758) Computational imaging; (170.1630) Coded aperture imaging;(100.5070) Phase retrieval; (090.1995) Digital holography.
References and links 1. B. Bhaduri, C. Edwards, H. Pham, R. Zhou, T. H. Nguyen, L. L. Goddard, and G. Popescu, “Diffraction phase microscopy: principles and applications in materials and life sciences,” Adv. Opt. Photon. 6, 57–119 (2014). 2. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996). 3. G. Nehmetallah and P. P. Banerjee, “Applications of digital and analog holography in three-dimensional imaging,” Adv. Opt. Photon. 4, 472–553 (2012). 4. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). 5. P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, “High-resolution scanning X-ray diffraction microscopy,” Science 321, 379–382 (2008). 6. N. Streibl, “Phase imaging by the transport equation of intensity,” Opt. Commun. 49, 6–10 (1984). 7. R. Horisaki, Y. Ogura, M. Aino, and J. Tanida, “Single-shot phase imaging with a coded aperture,” Opt. Lett. 39, 6466–6469 (2014). 8. Y. Rivenson, A. Stern, and B. Javidi, “Compressive Fresnel holography,” J. Disp. Technol. 6, 506–509 (2010). 9. D. L. Donoho, “Compressed sensing,” IEEE Trans. Info. Theory 52, 1289–1306 (2006). 10. D. J. Brady, K. Choi, D. L. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Express 17, 13040– 13049 (2009). 11. R. Horisaki, J. Tanida, A. Stern, and B. Javidi, “Multidimensional imaging using compressive Fresnel holography,” Opt. Lett. 37, 2013–2015 (2012). 12. M. Dierolf, A. Menzel, P. Thibault, P. Schneider, C. M. Kewish, R. Wepf, O. Bunk, and F. Pfeiffer, “Ptychographic X-ray computed tomography at the nanoscale,” Nature 467, 436–439 (2010). 13. D. J. Batey, D. Claus, and J. M. Rodenburg, “Information multiplexing in ptychography,” Ultramicroscopy 138, 13–21 (2014).
#232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9696
14. S. Dong, R. Shiradkar, P. Nanda, and G. Zheng, “Spectral multiplexing and coherent-state decomposition in Fourier ptychographic imaging,” Biomed. Opt. Express 5, 1757–1767 (2014). 15. T. M. Godden, R. Suman, M. J. Humphry, J. M. Rodenburg, and A. M. Maiden, “Ptychographic microscope for three-dimensional imaging,” Opt. Express 22, 12513–12523 (2014). 16. J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Trans. Image Proc. 16, 2992–3004 (2007). 17. A. M. Maiden, M. J. Humphry, and J. M. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29, 1606–1614 (2012). 18. L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica 2, 104–111 (2015). 19. “The USC-SIPI Image Database,” http://sipi.usc.edu/database/. 20. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D 60, 259–268 (1992). 21. Y. Rivenson and A. Stern, “Conditions for practicing compressive Fresnel holography,” Opt. Lett. 36, 3365–3367 (2011). 22. R. H. T. Bates, “Uniqueness of solutions to two-dimensional Fourier phase problems for localized and positive images,” Comput. Vision, Graph. Image Process. 25, 205–217 (1984). 23. N. Chacko, M. Liebling, and T. Blu, “Discretization of continuous convolution operators for accurate modeling of wave propagation in digital holography,” J. Opt. Soc. Am. A 30, 2012–2020 (2013).
1.
Introduction
Phase imaging is a powerful tool, especially for life science because many biological specimens have low light absorption and are difficult to detect using only the amplitude of propagating light coming from them [1]. A typical approach for phase imaging is an interferometric method such as holography [2, 3]. In general, this approach requires a complex setup for reference light or multiple shots for phase-shifting to quantitatively determine the interference. Another established approach for phase imaging without using light interference is coherent diffraction imaging (CDI). In CDI, a complex field of an object is reconstructed from a single intensity image of a diffraction pattern from the object with a phase retrieval (PR) algorithm [4]. To solve the reconstruction problem, spatial assumptions about the object, called support, are needed to constrain the solution. However, such assumptions severely restrict the size of the object. Ptychography and the transport-of-intensity equation (TIE) are famous multi-shot CDI methods for overcoming the object size limitation [5, 6]. We have recently proposed a single-shot CDI method for observing a large complex object. This method, referred to as single-shot phase imaging with a coded aperture (SPICA), has been used to demonstrate acquisition of two-dimensional complex objects [7]. SPICA connects synthetic aperture-based holography and CDI by means of a coded aperture (CA), which is a randomly arranged pinhole array that sieves the propagating field from the object, as shown in Fig. 1. The holographic process uses compressive Fresnel holography (CFH). In CFH, a complex object is reconstructed from randomly-sampled sparse holographic data by a numerical algorithm using the sparsity of the object based on compressive sensing (CS) [8]. CS is a recently-emerged sampling technique for observing object data with fewer measurements compared with the conventional sampling theorem [9]. We extend SPICA to single-shot acquisition of multidimensional complex objects in this paper. Several imaging methods for acquiring objects with more than two dimensions via holography or ptychography have been proposed [10–15]. Here we apply them to multidimensional SPICA. We employ an optical element, which is called coding optics in this paper, to modulate the sieved multidimensional complex field after the CA, as shown in Fig. 1. This allows us to realize single-shot multidimensional complex object acquisition with compact hardware that does not need reference light. A generalized framework for multidimensional SPICA is derived in Section 2, and its reconstruction process is shown in Section 3. Then, as an example, we focus on depth-variant
#232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9697
M chr o n o om e
Coding optics Multidimensional compressive Coherent diffraction imaging Fresnel holography with optical encoding
Fig. 1. Schematic diagram of multidimensional SPICA.
multispectral object acquisition based on the framework in Section 4, where its imaging performance is also demonstrated. A theoretical analysis of the proposed method is provided and numerically verified in Section 5. We summarize this paper in Section 6. 2.
Generalized system model
We show a generalized imaging process for multidimensional SPICA with matrix operators that are depicted in Fig. 1. First, we treat the matrix operators for describing the imaging process in an abstract form, and then provide their specific forms in the following part of this section and Section 4. In the process description, we refer to x-y planes in the multidimensional object as a channel. If individual channels along a dimensional axis have coherence between them, the dimension is called a coherent dimension, such as depth, when the depths are included within a coherence length. On the other hand, when individual channels along a dimensional axis do not have coherence between them, the dimension is defined as an incoherent dimension. For example, a spectrum is an incoherent dimension when each spectrum arises from a different illumination source. In the following discussion, we show a system model simply with one coherent dimension and one incoherent dimension. This model can be easily extended to one with multiple coherent and/or multiple incoherent dimensions by rearranging the dimensions. The y, t and v-dimensions that appear in Fig. 1 are omitted for simplicity. The system model of multidimensional SPICA is composed of two processes, CFH and CDI, which are shown in Fig. 1. The CFH process is written as a = M I coh P f .
(1)
First, a multidimensional complex field f ∈ C(Nx ×Nc ×Nc )×1 of the object illuminated by cocoh inc coh inc herent light propagates with a Fresnel diffraction operator P ∈ C(Nx ×Nc ×Nc )×(N x ×Nc ×Nc ) to the CA. Next, the propagating field is coherently multiplexed with an integration operainc coh inc tor I coh ∈ R(Nx ×Nc )×(N x ×Nc ×Nc ) along the coherent dimensional axis. Lastly, the multiplexed inc inc field is sieved by the CA with a masking operator M ∈ R(Nx ×Nc )×(N x ×Nc ) , as the sieved cominc plex field a ∈ C(Nx ×Nc )×1 just behind the CA. Here, N x , Nccoh , and Ncinc are the numbers of elements along the x-axis, and the coherent and incoherent dimensional axes, respectively. The CDI process is written as coh
E a |2 . g = I inc |E
#232104 - $15.00 USD © 2015 OSA
inc
(2)
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9698
First, channels of the sieved complex field a just behind the CA are modulated with an eninc inc coding operator E ∈ C(Nx ×Nc )×(N x ×Nc ) by the coding optics. The modulated field is then incoherently multiplexed and finally captured by a monochrome image sensor with an integrainc tion operator I inc ∈ RNx ×(Nx ×Nc ) along the incoherent dimensional axis, as the single intensity N ×1 image g ∈ R x . The encoding operator E can be arbitrarily chosen to modulate target dimensions, as shown in Section 4, but it needs to have phase dependencies for the PR process described in Section 3. In specific forms, Eq. (1), which shows the CFH process, is rewritten as ⎤ ⎡ ⎤ ⎤ ⎡ ⎡ 0 ··· 0 ⎥⎥⎥ 0 ⎥⎥ ⎢⎢⎢I coh 0 ··· 0 ⎥⎥⎥ ⎢⎢⎢P 1,1 ⎢⎢⎢ M 0 0 · · · 0 ⎥ ⎥ ⎥⎥⎥ ⎢ ⎢ ⎢⎢⎢ .. .. ⎥⎥⎥⎥ ⎢⎢⎢ .. ⎥⎥⎥ ⎢⎢⎢ ⎥⎥⎥ coh ⎢⎢⎢ 0 M ⎥ ⎢ ⎢ ⎥ 0 . 0 P ⎥ ⎢ ⎢ 0 . 0 . 0 I 2,1 ⎥ 0 ⎥⎥⎥ × ⎢⎢⎢ ⎥⎥⎥ f . (3) ⎥⎥⎥ × ⎢⎢⎢ 0 a = ⎢⎢⎢⎢ . ⎥ ⎥⎥⎥ ⎢ ⎢ ⎥⎥⎥ ⎢⎢ .. .. ⎥⎥⎥ ⎢⎢⎢ .. .. .. ⎢⎢⎢ . ⎥⎥⎥ ⎢ ⎥ . . . 0 0 0 0 ⎥⎥ ⎢⎢⎢ . 0 0 ⎥⎥⎥ ⎢⎢⎢ . ⎢⎣⎢ . ⎦ ⎣ ⎦ ⎣ ⎦⎥ coh 0 ··· 0 P Nccoh ,Ncinc 0 ··· 0 M0 0 ··· 0 I0 The masking matrix M in Eq. (1) is composed of a diagonal sub-matrix M 0 ∈ RNx ×Nx , where the coh the coherent diagonal elements correspond to the mask pattern. The integration matrix I for coh coh N ×(N x x ×Nc ) , where I ∈ dimension in Eq. (1) is composed of a sub-matrix I 0 = I 0 · · · I 0 ∈ R 0 RNx ×Nx is an identity matrix. The propagation matrix P in Eq. (1) is composed of Toeplitz sub-matrices PC coh ,C inc ∈ CNx ×Nx , which expresses the convolution of the Fresnel point spread function (PSF) at the C coh -th coherent and C inc -th incoherent channel. In the same manner, Eq. (2), which expresses the CDI process, is rewritten as
g = I0
I0
···
⎡
⎢⎢⎢ E 1
⎢⎢
⎢⎢⎢⎢ 0 I 0 ×
⎢⎢⎢⎢ .
⎢⎢⎢ .
⎢⎢⎢⎣ .
0
0
···
E2
0 .. . 0
0 ···
⎤
2 ⎥⎥⎥
⎥⎥⎥
⎥⎥⎥
⎥⎥⎥ a . ⎥⎥
0 ⎥⎥⎥⎥⎥
⎦
E N inc
0 .. .
(4)
c
The integration matrix I inc for the incoherent dimension in Eq. (2) is composed of the identity matrix I 0 . The encoding matrix E in Eq. (2) consists of sub-matrices E C inc ∈ CNx ×Nx , which causes each C inc -th channel of the sieved complex field a to be differently modulated by the coding optics. 3.
Reconstruction process
When the object is two-dimensional (Nccoh = 1, Ncinc = 1), the inversion of Eq. (2) is solved by the input-output PR (IOPR) algorithm with a CA-based support constraint [4, 7]. The modified CA-based IOPR for multidimensional SPICA is as follows. In the modified IOPR, the complex
T T inc = E a ∈ C(Nx ×Nc )×1 on the sensor plane before the incoherent field h = (hh1 )T · · · h N inc c
integration is also estimated. Here, hC inc ∈ CNx ×1 is the complex field of the C inc -th incoherent channel on the sensor plane, and the superscript T is the transpose of a matrix. 1. The initial complex field hˇ on the sensor plane is generated as a random complex field. 2. The new estimation hˆ of the complex field on the sensor plane is calculated as ⎛
⎞1
⎜⎜⎜
hˇ inc (U)
2 g (U) ⎟⎟⎟ 2
C ⎟⎟⎟ exp i arg hˇ inc (U) , hˆ C inc (U) = ⎜⎜⎜⎜⎝
⎟ C 2
ˇ
⎠ C inc h C inc (U)
#232104 - $15.00 USD © 2015 OSA
(5)
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9699
where w (J) shows the J-th element of the vector w , and arg(w) is the argument of a complex number w. 3. The complex field hˆ on the sensor plane is back-propagated to the CA plane as aˇ = E −1 hˆ ∈ inc C(Nx ×Nc )×1 . 4. The new estimation aˆ of the sieved complex field just behind the CA is calculated with the CA-based support constraint as aˆ = M aˇ . 5. The sieved complex field aˆ is forward-propagated to the sensor plane as hˇ = E aˆ . Steps 2 to 5 are iterated until aˆ converges. Since the inversion of Eq. (1) is ill-posed, we use a sparsity-based CS algorithm called twostep iterative shrinkage/thresholding (TwIST) for the inversion [16]. TwIST solves the following optimization problem iteratively: ˆf = arg min ||aˆ − M I coh P f || + τR( f ), 2
(6)
f
where || • ||2 is the 2 norm, τ is a regularization parameter, and R(•) is a regularizer. In the system model and the reconstruction of the CFH process, we used the Born approximation (single scattering model) [2, 10]. This approximation degrades the performance of the reconstruction when the object is depth-variant and spatially dense. This issue can be solved by employing some state-of-the-art algorithms used in ptychography [15, 17, 18]. 4.
Depth-variant multispectral object acquisition
Here, as a special case of Section 2, we show single-shot acquisition of a depth-variant multispectral (four-dimensional) complex object F (x, y, z, λ), where x, y, and z are the three spatial dimensions, as shown in Fig. 1, and λ represents the spectral dimension. The depth z and the spectrum λ are the coherent and incoherent dimensions, respectively. To modulate the two dimensions z and λ, as the coding optics we choose free space , described with the encoding operator E in Eq. (2). Free space results in Fresnel propagation, and this propagation process is dependent on the depth z, the wavelength λ, and also the phase. This imaging modality is readily applicable to other dimensions, such as time and polarization. The setup for SPICA-based depth-variant multispectral object acquisition is composed of a CA, a monochrome image sensor, and coherent illumination. The mathematical model of this system is written as A(s, λ) = M(s) P(s − x, z, λ)F (x, z, λ)dxdz, (7)
2 (8) G(u) =
P(u − s, zA , λ)A(s, λ)ds
dλ, where P(•, z, λ) is the PSF of the Fresnel propagating light at distance z and wavelength λ, M is the mask pattern of the CA, A is the sieved complex field just behind the CA, and G is the single intensity image captured by the monochrome image sensor. s and u are the spatial dimensions on the CA and sensor planes, and zA is the distance between the CA and sensor planes, as shown in Fig. 1. Equation (7) expresses the CFH process shown in Eq. (1). The depth-variant multispectral complex field F of the object illuminated by polychromatic coherent light propagates with the Fresnel diffraction kernel P to the CA, the propagating field is coherently multiplexed along the z-axis, and the multiplexed field is multiplied with the mask pattern M of the CA, resulting in the sieved complex field A. Equation (8) corresponds to the CDI process shown in Eq. (2). #232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9700
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
Fig. 2. Simulation of multidimensional SPICA. The (a) amplitude and (b) phase of a color object at the distance z1 . The (c) amplitude and (d) phase of the color object at the distance z2 . (e) A magnified image (pixel count, 30 × 30) of the coded aperture. (f) The captured intensity image. The (g) amplitude and (h) phase of the reconstruction at the distance z1 . The (i) amplitude and (j) phase of the reconstruction at the distance z2 . Phases are normalized in the interval [−π, π].
The sieved complex field A just behind the CA propagates with the Fresnel diffraction kernel P to the monochrome image sensor. The propagating field is incoherently multiplexed along the λ-axis onto the monochrome image sensor, generating the captured single intensity image G. As indicated in the previous paragraph, the imaging process described by Eqs. (7) and (8) is a special case of the generalized system model of multidimensional SPICA described with Eqs. (1) and (2). Here Nccoh = Nz is the number of elements along the z-axis, and Ncinc = Nλ is the number of elements along the λ-axis. Also the Toeplitz sub-matrix PC coh ,C inc = P Z,Λ expresses the convolution of the Fresnel PSF at the Z-th depth and Λ-th spectral channel, in Eq. (3), and the encoding sub-matrix E C inc = P zΛA also expresses the convolution of the Fresnel PSF at distance zA and the Λ-th spectral channel, in Eq. (4), respectively. We verified SPICA-based single-shot depth-variant multispectral object acquisition numerically, as shown in Fig. 2. In the simulation, the depth-variant multispectral complex object f shown in Figs. 2(a)–2(d) was composed of two depths (Nz = 2) at distances 30 mm (z1 ) and 60 mm (z2 ) as the coherent dimension from the CA and three spectrums (Nλ = 3), namely, 460 nm (λ1 ), 530 nm (λ2 ), and 630 nm (λ3 ) as the incoherent dimension. The channels consisted of natural images [19]. The pixel count N x 2 was 5122 , and the pixel pitch δ x was 5 μm. A CA whose magnified mask pattern M 0 is shown in Figs. 2(e) was located 3 mm (zA ) from the sensor plane. The number NP of the randomly arranged pinholes, shown as white pixels in the figure, was 21845, which was 8.3 % of the total pixel count. The captured intensity image g is shown in Fig. 2(f), where Gaussian random noise with a signal-to-noise ratio (SNR) of 30 dB was added. The reconstruction result ˆf after phase bias matching is shown in Figs. 2(g)–2(j). We used the two-dimensional total variation as the regularizer R in Eq. (6) for TwIST [20]. The complex object field was reconstructed successfully, and the peak SNR (PSNR) was 29.5 dB. The effects of the measurement noise level and the alignment error of the pinholes in the CA on the reconstruction fidelity are shown in Fig. 3. The former effect is shown in Fig. 3(a), in which the relationship between the measurement SNR and the the reconstruction PSNR are #232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9701
35 Reconstruction PSNR (dB)
Reconstruction PSNR (dB)
35 30 25 20 15 10 5 0 0
10 20 30 40 Measurement SNR (dB)
50
30 25 20 15 10 5 0 0
(a)
0.1
0.2 0.3 Support error
0.4
0.5
(b)
Fig. 3. Effects for the reconstruction fidelity by (a) the measurement noise level and (b) the error of the CA-based support constraint.
plotted. The latter effect is shown in Fig. 3(b). In this figure, to emulate the alignment error, the locations of multiple pinholes in the CA-based support constraint in the reconstruction process are randomly shifted by one pixel from those on the CA in the simulated imaging process. The ratio of the number of shifted pinholes to the total number of pinholes is defined as the support error. The measurement SNR was fixed to 30 dB in Fig. 3(b). The simulated setup used for obtaining Fig. 3 is the same as that used for the previous simulation in Fig. 2. Those plots show averaged reconstruction PSNRs of ten different mask patterns. The reconstruction PSNRs rise at a measurement SNR of about 20 dB and reach a lower bound at around the support error of 0.25 in this case. 5.
Theoretical analysis
The theoretical conditions for multidimensional SPICA are derived from the CFH process in Eq. (1) and the CDI process in Eq. (2) [7]. The following discussions assume a depth-variant multispectral complex object as in the previous section. Those discussions are extendable to other dimensions, besides the system model described in Section 2. When the object distance z from the CA satisfies z ≥ zmin =
Nx δx 2 , λ
(9)
CFH is approximated as an ideal three-dimensional Fourier encoder [10]. Then the number NCFH of the required sampling points in each spectral channel in CFH is calculated, based on a worst-case analysis, to be
σK log N x 2 Nz , (10) NCFH ≥ Nλ where σ is a constant depending on the acceptable reconstruction fidelity, and K is the sparsity in the regularizer domain for all depth channels and spectral channels [21]. In the SPICA case, NCFH = NP , which is the number of the pinholes in the CA. In the previous simulation, zmin in Eq. (9) with the minimum value of the wavelength λ was 28 mm, and thus the condition was satisfied. #232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9702
In the CDI process, oversampling is required to find unique autocorrelation of a complex object [4]. In multidimensional SPICA, the oversampling factor may be Nλ -times larger for PR in each spectral channel compared with that of two-dimensional CDI. The distance zA between the CA and sensor planes needs to be √ ρ2D Nλ δ x 2 , (11) = zA ≥ zmin A λ where ρ2D is the oversampling factor of PR for a two-dimensional object, because a single pixel on the CA contributes more than ρ2D Nλ pixels on the image sensor in this condition. The upper bound of the number NP of pinholes on the CA is calculated as NP ≤ NPmax =
Nx 2 . ρ2D Nλ
(12)
The two-dimensional oversampling factor is known to be ρ2D ≈ 4 [22]. In the previous simumax in lation, zmin A in Eq. (11) with the minimum value of the wavelength λ was 0.2 mm and NP Eq. (12) was 21845, and thus the conditions were satisfied. The upper bound of multidimensional SPICA is calculated with Eqs. (10) and (12) to be σK 1
≥ 2 = σq, 2 Nx ρ2D log N x Nz
(13)
where q = K/N x 2 is defined as the relative sparsity. This condition indicates that the upper bound of the relative sparsity q is roughly independent for the numbers Nz and Nλ of depth channels and spectral channels when the pixel count Nx 2 of the image sensor is sufficiently large. The lateral and axial resolutions of multidimensional SPICA are equivalent to those in conand δmin are calculated ventional holographic imaging [10]. The achievable resolutions δmin x z as λz , Nx δx 4λz2 δmin = 2 2. z Nx δx
δmin x =
(14) (15)
The spectral resolution of multidimensional SPICA is estimated with the Airy disk on the sensor plane using a pinhole of diameter δ x on the CA. The diameter d of the first dark annulus of the Airy disk is calculated to be d = 2.44λzA /δ x [2], and ∂d/∂λ = 2.44zA /δ x . To change the diameter d by more than one pixel pitch δ x on the image sensor, the achievable spectral resolution δmin λ is calculated to be δx 2 . (16) δmin λ = 2.44zA The theoretical conditions discussed above were verified numerically as shown in Fig. 4. The plot shows the relationship between the required relative sparsity q and the numbers Nz and Nλ of depth channels and spectral channels. The object was composed of point sources randomly distributed four-dimensionally and having the same intensity and different random phases. The sparsity K is the number of sources in this case. The pixel count Nx 2 was 642 , the pixel pitch δ x was 5 μm, and the minimum wavelength λ was 530 nm. The nearest plane was located at zmin (= 3 mm) in Eq. (9) with the minimum wavelength. The distance zA of the CA was set to zmin A (a function of Nλ ) in Eq. (11) with the minimum wavelength. The number NP #232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9703
0.08 0.06 0.04 0.02 0 1
2
3
4 1
2
3
4
Fig. 4. Relationship between the relative sparsity q and the numbers Nz and Nλ of depth channels and spectral channels. The numerically and analytically calculated relative sparsities q are shown with the blue dots and the red grid, respectively.
of pinholes on the CA was chosen to be NPmax (a function of Nλ ) in Eq (12). The axial pitch in Eq. (15) with the maximum distance (a function of Nz ) and the maximum was set to be δmin z wavelength (a function of Nλ ). The spectral pitch was set to δmin λ (a function of Nλ ) in Eq. (16). In Fig. 4, each of the numerically calculated relative sparsities q, which assume a PSNR of 30 dB as the acceptable reconstruction fidelity, are shown with the blue dots. They are averages obtained using ten combinations of different mask patterns of the CA and distributions of the point sources. The 1 norm was chosen as the regularizer R in Eq. (6) for TwIST. The analytically calculated relative sparsities q shown with the red grid in Fig. 4 were from Eq. (13) with σ = 0.6, which was chosen experimentally. The numerical relative sparsity q approaches the analytical one, which is approximately independent of the numbers Nz and Nλ of depth channels and spectral channels in this case, as N x and Nλ increase, although there is a small difference between them based on the worst-case analysis. 6.
Conclusion
In summary, we extended SPICA to single-shot acquisition of multidimensional complex objects. The generalized system model and its reconstruction process were presented. As an example, assuming a depth-variant multispectral object, the system model and the theoretical conditions were derived and numerically demonstrated. In the multidimensional SPICA method, both the depth and spectrum dimensions are encoded by free space propagation, and the modulated depth and spectral dimensions are coherently and incoherently multiplexed, respectively. In the numerical demonstrations, we showed the effect of the error between the CA-based support constraint and the actual mask pattern in the reconstruction process. Other modeling errors, such as discretization [23], should also be investigated. The proposed imaging modality can be readily applied to other optical dimensions, such as time and polarization. Our approach consists of only a CA, a monochrome image sensor, and a coherent light source, and does not require any reference light. Thus, a compact and highly functional imaging system can be realized. In addition, the technique is applicable to incoherent holography and sensing in various spectral regions, such as X-ray and electron imaging.
#232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9704
Abstract: We propose a generalized framework for quantitatively acquiring multidimensional complex objects based on single-shot phase imaging with a coded aperture (SPICA). In multidimensional SPICA, a propagating field from a multidimensional complex object is sieved by a coded aperture, the sieved field is modulated by an optical element, which is called coding optics, and then the resultant field is captured by a monochrome image sensor. The original complex field is reconstructed from the single captured intensity image by a phase retrieval algorithm with a support constraint of the coded aperture and a sparsity-based reconstruction algorithm based on compressive sensing. We also present theoretical conditions for the proposed method. As a demonstration, we numerically verified an application of this generalized framework for single-shot acquisition of depth-variant multispectral objects. © 2015 Optical Society of America OCIS codes: (110.1758) Computational imaging; (170.1630) Coded aperture imaging;(100.5070) Phase retrieval; (090.1995) Digital holography.
References and links 1. B. Bhaduri, C. Edwards, H. Pham, R. Zhou, T. H. Nguyen, L. L. Goddard, and G. Popescu, “Diffraction phase microscopy: principles and applications in materials and life sciences,” Adv. Opt. Photon. 6, 57–119 (2014). 2. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996). 3. G. Nehmetallah and P. P. Banerjee, “Applications of digital and analog holography in three-dimensional imaging,” Adv. Opt. Photon. 4, 472–553 (2012). 4. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). 5. P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, “High-resolution scanning X-ray diffraction microscopy,” Science 321, 379–382 (2008). 6. N. Streibl, “Phase imaging by the transport equation of intensity,” Opt. Commun. 49, 6–10 (1984). 7. R. Horisaki, Y. Ogura, M. Aino, and J. Tanida, “Single-shot phase imaging with a coded aperture,” Opt. Lett. 39, 6466–6469 (2014). 8. Y. Rivenson, A. Stern, and B. Javidi, “Compressive Fresnel holography,” J. Disp. Technol. 6, 506–509 (2010). 9. D. L. Donoho, “Compressed sensing,” IEEE Trans. Info. Theory 52, 1289–1306 (2006). 10. D. J. Brady, K. Choi, D. L. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Express 17, 13040– 13049 (2009). 11. R. Horisaki, J. Tanida, A. Stern, and B. Javidi, “Multidimensional imaging using compressive Fresnel holography,” Opt. Lett. 37, 2013–2015 (2012). 12. M. Dierolf, A. Menzel, P. Thibault, P. Schneider, C. M. Kewish, R. Wepf, O. Bunk, and F. Pfeiffer, “Ptychographic X-ray computed tomography at the nanoscale,” Nature 467, 436–439 (2010). 13. D. J. Batey, D. Claus, and J. M. Rodenburg, “Information multiplexing in ptychography,” Ultramicroscopy 138, 13–21 (2014).
#232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9696
14. S. Dong, R. Shiradkar, P. Nanda, and G. Zheng, “Spectral multiplexing and coherent-state decomposition in Fourier ptychographic imaging,” Biomed. Opt. Express 5, 1757–1767 (2014). 15. T. M. Godden, R. Suman, M. J. Humphry, J. M. Rodenburg, and A. M. Maiden, “Ptychographic microscope for three-dimensional imaging,” Opt. Express 22, 12513–12523 (2014). 16. J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Trans. Image Proc. 16, 2992–3004 (2007). 17. A. M. Maiden, M. J. Humphry, and J. M. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29, 1606–1614 (2012). 18. L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica 2, 104–111 (2015). 19. “The USC-SIPI Image Database,” http://sipi.usc.edu/database/. 20. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D 60, 259–268 (1992). 21. Y. Rivenson and A. Stern, “Conditions for practicing compressive Fresnel holography,” Opt. Lett. 36, 3365–3367 (2011). 22. R. H. T. Bates, “Uniqueness of solutions to two-dimensional Fourier phase problems for localized and positive images,” Comput. Vision, Graph. Image Process. 25, 205–217 (1984). 23. N. Chacko, M. Liebling, and T. Blu, “Discretization of continuous convolution operators for accurate modeling of wave propagation in digital holography,” J. Opt. Soc. Am. A 30, 2012–2020 (2013).
1.
Introduction
Phase imaging is a powerful tool, especially for life science because many biological specimens have low light absorption and are difficult to detect using only the amplitude of propagating light coming from them [1]. A typical approach for phase imaging is an interferometric method such as holography [2, 3]. In general, this approach requires a complex setup for reference light or multiple shots for phase-shifting to quantitatively determine the interference. Another established approach for phase imaging without using light interference is coherent diffraction imaging (CDI). In CDI, a complex field of an object is reconstructed from a single intensity image of a diffraction pattern from the object with a phase retrieval (PR) algorithm [4]. To solve the reconstruction problem, spatial assumptions about the object, called support, are needed to constrain the solution. However, such assumptions severely restrict the size of the object. Ptychography and the transport-of-intensity equation (TIE) are famous multi-shot CDI methods for overcoming the object size limitation [5, 6]. We have recently proposed a single-shot CDI method for observing a large complex object. This method, referred to as single-shot phase imaging with a coded aperture (SPICA), has been used to demonstrate acquisition of two-dimensional complex objects [7]. SPICA connects synthetic aperture-based holography and CDI by means of a coded aperture (CA), which is a randomly arranged pinhole array that sieves the propagating field from the object, as shown in Fig. 1. The holographic process uses compressive Fresnel holography (CFH). In CFH, a complex object is reconstructed from randomly-sampled sparse holographic data by a numerical algorithm using the sparsity of the object based on compressive sensing (CS) [8]. CS is a recently-emerged sampling technique for observing object data with fewer measurements compared with the conventional sampling theorem [9]. We extend SPICA to single-shot acquisition of multidimensional complex objects in this paper. Several imaging methods for acquiring objects with more than two dimensions via holography or ptychography have been proposed [10–15]. Here we apply them to multidimensional SPICA. We employ an optical element, which is called coding optics in this paper, to modulate the sieved multidimensional complex field after the CA, as shown in Fig. 1. This allows us to realize single-shot multidimensional complex object acquisition with compact hardware that does not need reference light. A generalized framework for multidimensional SPICA is derived in Section 2, and its reconstruction process is shown in Section 3. Then, as an example, we focus on depth-variant
#232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9697
M chr o n o om e
Coding optics Multidimensional compressive Coherent diffraction imaging Fresnel holography with optical encoding
Fig. 1. Schematic diagram of multidimensional SPICA.
multispectral object acquisition based on the framework in Section 4, where its imaging performance is also demonstrated. A theoretical analysis of the proposed method is provided and numerically verified in Section 5. We summarize this paper in Section 6. 2.
Generalized system model
We show a generalized imaging process for multidimensional SPICA with matrix operators that are depicted in Fig. 1. First, we treat the matrix operators for describing the imaging process in an abstract form, and then provide their specific forms in the following part of this section and Section 4. In the process description, we refer to x-y planes in the multidimensional object as a channel. If individual channels along a dimensional axis have coherence between them, the dimension is called a coherent dimension, such as depth, when the depths are included within a coherence length. On the other hand, when individual channels along a dimensional axis do not have coherence between them, the dimension is defined as an incoherent dimension. For example, a spectrum is an incoherent dimension when each spectrum arises from a different illumination source. In the following discussion, we show a system model simply with one coherent dimension and one incoherent dimension. This model can be easily extended to one with multiple coherent and/or multiple incoherent dimensions by rearranging the dimensions. The y, t and v-dimensions that appear in Fig. 1 are omitted for simplicity. The system model of multidimensional SPICA is composed of two processes, CFH and CDI, which are shown in Fig. 1. The CFH process is written as a = M I coh P f .
(1)
First, a multidimensional complex field f ∈ C(Nx ×Nc ×Nc )×1 of the object illuminated by cocoh inc coh inc herent light propagates with a Fresnel diffraction operator P ∈ C(Nx ×Nc ×Nc )×(N x ×Nc ×Nc ) to the CA. Next, the propagating field is coherently multiplexed with an integration operainc coh inc tor I coh ∈ R(Nx ×Nc )×(N x ×Nc ×Nc ) along the coherent dimensional axis. Lastly, the multiplexed inc inc field is sieved by the CA with a masking operator M ∈ R(Nx ×Nc )×(N x ×Nc ) , as the sieved cominc plex field a ∈ C(Nx ×Nc )×1 just behind the CA. Here, N x , Nccoh , and Ncinc are the numbers of elements along the x-axis, and the coherent and incoherent dimensional axes, respectively. The CDI process is written as coh
E a |2 . g = I inc |E
#232104 - $15.00 USD © 2015 OSA
inc
(2)
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9698
First, channels of the sieved complex field a just behind the CA are modulated with an eninc inc coding operator E ∈ C(Nx ×Nc )×(N x ×Nc ) by the coding optics. The modulated field is then incoherently multiplexed and finally captured by a monochrome image sensor with an integrainc tion operator I inc ∈ RNx ×(Nx ×Nc ) along the incoherent dimensional axis, as the single intensity N ×1 image g ∈ R x . The encoding operator E can be arbitrarily chosen to modulate target dimensions, as shown in Section 4, but it needs to have phase dependencies for the PR process described in Section 3. In specific forms, Eq. (1), which shows the CFH process, is rewritten as ⎤ ⎡ ⎤ ⎤ ⎡ ⎡ 0 ··· 0 ⎥⎥⎥ 0 ⎥⎥ ⎢⎢⎢I coh 0 ··· 0 ⎥⎥⎥ ⎢⎢⎢P 1,1 ⎢⎢⎢ M 0 0 · · · 0 ⎥ ⎥ ⎥⎥⎥ ⎢ ⎢ ⎢⎢⎢ .. .. ⎥⎥⎥⎥ ⎢⎢⎢ .. ⎥⎥⎥ ⎢⎢⎢ ⎥⎥⎥ coh ⎢⎢⎢ 0 M ⎥ ⎢ ⎢ ⎥ 0 . 0 P ⎥ ⎢ ⎢ 0 . 0 . 0 I 2,1 ⎥ 0 ⎥⎥⎥ × ⎢⎢⎢ ⎥⎥⎥ f . (3) ⎥⎥⎥ × ⎢⎢⎢ 0 a = ⎢⎢⎢⎢ . ⎥ ⎥⎥⎥ ⎢ ⎢ ⎥⎥⎥ ⎢⎢ .. .. ⎥⎥⎥ ⎢⎢⎢ .. .. .. ⎢⎢⎢ . ⎥⎥⎥ ⎢ ⎥ . . . 0 0 0 0 ⎥⎥ ⎢⎢⎢ . 0 0 ⎥⎥⎥ ⎢⎢⎢ . ⎢⎣⎢ . ⎦ ⎣ ⎦ ⎣ ⎦⎥ coh 0 ··· 0 P Nccoh ,Ncinc 0 ··· 0 M0 0 ··· 0 I0 The masking matrix M in Eq. (1) is composed of a diagonal sub-matrix M 0 ∈ RNx ×Nx , where the coh the coherent diagonal elements correspond to the mask pattern. The integration matrix I for coh coh N ×(N x x ×Nc ) , where I ∈ dimension in Eq. (1) is composed of a sub-matrix I 0 = I 0 · · · I 0 ∈ R 0 RNx ×Nx is an identity matrix. The propagation matrix P in Eq. (1) is composed of Toeplitz sub-matrices PC coh ,C inc ∈ CNx ×Nx , which expresses the convolution of the Fresnel point spread function (PSF) at the C coh -th coherent and C inc -th incoherent channel. In the same manner, Eq. (2), which expresses the CDI process, is rewritten as
g = I0
I0
···
⎡
⎢⎢⎢ E 1
⎢⎢
⎢⎢⎢⎢ 0 I 0 ×
⎢⎢⎢⎢ .
⎢⎢⎢ .
⎢⎢⎢⎣ .
0
0
···
E2
0 .. . 0
0 ···
⎤
2 ⎥⎥⎥
⎥⎥⎥
⎥⎥⎥
⎥⎥⎥ a . ⎥⎥
0 ⎥⎥⎥⎥⎥
⎦
E N inc
0 .. .
(4)
c
The integration matrix I inc for the incoherent dimension in Eq. (2) is composed of the identity matrix I 0 . The encoding matrix E in Eq. (2) consists of sub-matrices E C inc ∈ CNx ×Nx , which causes each C inc -th channel of the sieved complex field a to be differently modulated by the coding optics. 3.
Reconstruction process
When the object is two-dimensional (Nccoh = 1, Ncinc = 1), the inversion of Eq. (2) is solved by the input-output PR (IOPR) algorithm with a CA-based support constraint [4, 7]. The modified CA-based IOPR for multidimensional SPICA is as follows. In the modified IOPR, the complex
T T inc = E a ∈ C(Nx ×Nc )×1 on the sensor plane before the incoherent field h = (hh1 )T · · · h N inc c
integration is also estimated. Here, hC inc ∈ CNx ×1 is the complex field of the C inc -th incoherent channel on the sensor plane, and the superscript T is the transpose of a matrix. 1. The initial complex field hˇ on the sensor plane is generated as a random complex field. 2. The new estimation hˆ of the complex field on the sensor plane is calculated as ⎛
⎞1
⎜⎜⎜
hˇ inc (U)
2 g (U) ⎟⎟⎟ 2
C ⎟⎟⎟ exp i arg hˇ inc (U) , hˆ C inc (U) = ⎜⎜⎜⎜⎝
⎟ C 2
ˇ
⎠ C inc h C inc (U)
#232104 - $15.00 USD © 2015 OSA
(5)
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9699
where w (J) shows the J-th element of the vector w , and arg(w) is the argument of a complex number w. 3. The complex field hˆ on the sensor plane is back-propagated to the CA plane as aˇ = E −1 hˆ ∈ inc C(Nx ×Nc )×1 . 4. The new estimation aˆ of the sieved complex field just behind the CA is calculated with the CA-based support constraint as aˆ = M aˇ . 5. The sieved complex field aˆ is forward-propagated to the sensor plane as hˇ = E aˆ . Steps 2 to 5 are iterated until aˆ converges. Since the inversion of Eq. (1) is ill-posed, we use a sparsity-based CS algorithm called twostep iterative shrinkage/thresholding (TwIST) for the inversion [16]. TwIST solves the following optimization problem iteratively: ˆf = arg min ||aˆ − M I coh P f || + τR( f ), 2
(6)
f
where || • ||2 is the 2 norm, τ is a regularization parameter, and R(•) is a regularizer. In the system model and the reconstruction of the CFH process, we used the Born approximation (single scattering model) [2, 10]. This approximation degrades the performance of the reconstruction when the object is depth-variant and spatially dense. This issue can be solved by employing some state-of-the-art algorithms used in ptychography [15, 17, 18]. 4.
Depth-variant multispectral object acquisition
Here, as a special case of Section 2, we show single-shot acquisition of a depth-variant multispectral (four-dimensional) complex object F (x, y, z, λ), where x, y, and z are the three spatial dimensions, as shown in Fig. 1, and λ represents the spectral dimension. The depth z and the spectrum λ are the coherent and incoherent dimensions, respectively. To modulate the two dimensions z and λ, as the coding optics we choose free space , described with the encoding operator E in Eq. (2). Free space results in Fresnel propagation, and this propagation process is dependent on the depth z, the wavelength λ, and also the phase. This imaging modality is readily applicable to other dimensions, such as time and polarization. The setup for SPICA-based depth-variant multispectral object acquisition is composed of a CA, a monochrome image sensor, and coherent illumination. The mathematical model of this system is written as A(s, λ) = M(s) P(s − x, z, λ)F (x, z, λ)dxdz, (7)
2 (8) G(u) =
P(u − s, zA , λ)A(s, λ)ds
dλ, where P(•, z, λ) is the PSF of the Fresnel propagating light at distance z and wavelength λ, M is the mask pattern of the CA, A is the sieved complex field just behind the CA, and G is the single intensity image captured by the monochrome image sensor. s and u are the spatial dimensions on the CA and sensor planes, and zA is the distance between the CA and sensor planes, as shown in Fig. 1. Equation (7) expresses the CFH process shown in Eq. (1). The depth-variant multispectral complex field F of the object illuminated by polychromatic coherent light propagates with the Fresnel diffraction kernel P to the CA, the propagating field is coherently multiplexed along the z-axis, and the multiplexed field is multiplied with the mask pattern M of the CA, resulting in the sieved complex field A. Equation (8) corresponds to the CDI process shown in Eq. (2). #232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9700
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
Fig. 2. Simulation of multidimensional SPICA. The (a) amplitude and (b) phase of a color object at the distance z1 . The (c) amplitude and (d) phase of the color object at the distance z2 . (e) A magnified image (pixel count, 30 × 30) of the coded aperture. (f) The captured intensity image. The (g) amplitude and (h) phase of the reconstruction at the distance z1 . The (i) amplitude and (j) phase of the reconstruction at the distance z2 . Phases are normalized in the interval [−π, π].
The sieved complex field A just behind the CA propagates with the Fresnel diffraction kernel P to the monochrome image sensor. The propagating field is incoherently multiplexed along the λ-axis onto the monochrome image sensor, generating the captured single intensity image G. As indicated in the previous paragraph, the imaging process described by Eqs. (7) and (8) is a special case of the generalized system model of multidimensional SPICA described with Eqs. (1) and (2). Here Nccoh = Nz is the number of elements along the z-axis, and Ncinc = Nλ is the number of elements along the λ-axis. Also the Toeplitz sub-matrix PC coh ,C inc = P Z,Λ expresses the convolution of the Fresnel PSF at the Z-th depth and Λ-th spectral channel, in Eq. (3), and the encoding sub-matrix E C inc = P zΛA also expresses the convolution of the Fresnel PSF at distance zA and the Λ-th spectral channel, in Eq. (4), respectively. We verified SPICA-based single-shot depth-variant multispectral object acquisition numerically, as shown in Fig. 2. In the simulation, the depth-variant multispectral complex object f shown in Figs. 2(a)–2(d) was composed of two depths (Nz = 2) at distances 30 mm (z1 ) and 60 mm (z2 ) as the coherent dimension from the CA and three spectrums (Nλ = 3), namely, 460 nm (λ1 ), 530 nm (λ2 ), and 630 nm (λ3 ) as the incoherent dimension. The channels consisted of natural images [19]. The pixel count N x 2 was 5122 , and the pixel pitch δ x was 5 μm. A CA whose magnified mask pattern M 0 is shown in Figs. 2(e) was located 3 mm (zA ) from the sensor plane. The number NP of the randomly arranged pinholes, shown as white pixels in the figure, was 21845, which was 8.3 % of the total pixel count. The captured intensity image g is shown in Fig. 2(f), where Gaussian random noise with a signal-to-noise ratio (SNR) of 30 dB was added. The reconstruction result ˆf after phase bias matching is shown in Figs. 2(g)–2(j). We used the two-dimensional total variation as the regularizer R in Eq. (6) for TwIST [20]. The complex object field was reconstructed successfully, and the peak SNR (PSNR) was 29.5 dB. The effects of the measurement noise level and the alignment error of the pinholes in the CA on the reconstruction fidelity are shown in Fig. 3. The former effect is shown in Fig. 3(a), in which the relationship between the measurement SNR and the the reconstruction PSNR are #232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9701
35 Reconstruction PSNR (dB)
Reconstruction PSNR (dB)
35 30 25 20 15 10 5 0 0
10 20 30 40 Measurement SNR (dB)
50
30 25 20 15 10 5 0 0
(a)
0.1
0.2 0.3 Support error
0.4
0.5
(b)
Fig. 3. Effects for the reconstruction fidelity by (a) the measurement noise level and (b) the error of the CA-based support constraint.
plotted. The latter effect is shown in Fig. 3(b). In this figure, to emulate the alignment error, the locations of multiple pinholes in the CA-based support constraint in the reconstruction process are randomly shifted by one pixel from those on the CA in the simulated imaging process. The ratio of the number of shifted pinholes to the total number of pinholes is defined as the support error. The measurement SNR was fixed to 30 dB in Fig. 3(b). The simulated setup used for obtaining Fig. 3 is the same as that used for the previous simulation in Fig. 2. Those plots show averaged reconstruction PSNRs of ten different mask patterns. The reconstruction PSNRs rise at a measurement SNR of about 20 dB and reach a lower bound at around the support error of 0.25 in this case. 5.
Theoretical analysis
The theoretical conditions for multidimensional SPICA are derived from the CFH process in Eq. (1) and the CDI process in Eq. (2) [7]. The following discussions assume a depth-variant multispectral complex object as in the previous section. Those discussions are extendable to other dimensions, besides the system model described in Section 2. When the object distance z from the CA satisfies z ≥ zmin =
Nx δx 2 , λ
(9)
CFH is approximated as an ideal three-dimensional Fourier encoder [10]. Then the number NCFH of the required sampling points in each spectral channel in CFH is calculated, based on a worst-case analysis, to be
σK log N x 2 Nz , (10) NCFH ≥ Nλ where σ is a constant depending on the acceptable reconstruction fidelity, and K is the sparsity in the regularizer domain for all depth channels and spectral channels [21]. In the SPICA case, NCFH = NP , which is the number of the pinholes in the CA. In the previous simulation, zmin in Eq. (9) with the minimum value of the wavelength λ was 28 mm, and thus the condition was satisfied. #232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9702
In the CDI process, oversampling is required to find unique autocorrelation of a complex object [4]. In multidimensional SPICA, the oversampling factor may be Nλ -times larger for PR in each spectral channel compared with that of two-dimensional CDI. The distance zA between the CA and sensor planes needs to be √ ρ2D Nλ δ x 2 , (11) = zA ≥ zmin A λ where ρ2D is the oversampling factor of PR for a two-dimensional object, because a single pixel on the CA contributes more than ρ2D Nλ pixels on the image sensor in this condition. The upper bound of the number NP of pinholes on the CA is calculated as NP ≤ NPmax =
Nx 2 . ρ2D Nλ
(12)
The two-dimensional oversampling factor is known to be ρ2D ≈ 4 [22]. In the previous simumax in lation, zmin A in Eq. (11) with the minimum value of the wavelength λ was 0.2 mm and NP Eq. (12) was 21845, and thus the conditions were satisfied. The upper bound of multidimensional SPICA is calculated with Eqs. (10) and (12) to be σK 1
≥ 2 = σq, 2 Nx ρ2D log N x Nz
(13)
where q = K/N x 2 is defined as the relative sparsity. This condition indicates that the upper bound of the relative sparsity q is roughly independent for the numbers Nz and Nλ of depth channels and spectral channels when the pixel count Nx 2 of the image sensor is sufficiently large. The lateral and axial resolutions of multidimensional SPICA are equivalent to those in conand δmin are calculated ventional holographic imaging [10]. The achievable resolutions δmin x z as λz , Nx δx 4λz2 δmin = 2 2. z Nx δx
δmin x =
(14) (15)
The spectral resolution of multidimensional SPICA is estimated with the Airy disk on the sensor plane using a pinhole of diameter δ x on the CA. The diameter d of the first dark annulus of the Airy disk is calculated to be d = 2.44λzA /δ x [2], and ∂d/∂λ = 2.44zA /δ x . To change the diameter d by more than one pixel pitch δ x on the image sensor, the achievable spectral resolution δmin λ is calculated to be δx 2 . (16) δmin λ = 2.44zA The theoretical conditions discussed above were verified numerically as shown in Fig. 4. The plot shows the relationship between the required relative sparsity q and the numbers Nz and Nλ of depth channels and spectral channels. The object was composed of point sources randomly distributed four-dimensionally and having the same intensity and different random phases. The sparsity K is the number of sources in this case. The pixel count Nx 2 was 642 , the pixel pitch δ x was 5 μm, and the minimum wavelength λ was 530 nm. The nearest plane was located at zmin (= 3 mm) in Eq. (9) with the minimum wavelength. The distance zA of the CA was set to zmin A (a function of Nλ ) in Eq. (11) with the minimum wavelength. The number NP #232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9703
0.08 0.06 0.04 0.02 0 1
2
3
4 1
2
3
4
Fig. 4. Relationship between the relative sparsity q and the numbers Nz and Nλ of depth channels and spectral channels. The numerically and analytically calculated relative sparsities q are shown with the blue dots and the red grid, respectively.
of pinholes on the CA was chosen to be NPmax (a function of Nλ ) in Eq (12). The axial pitch in Eq. (15) with the maximum distance (a function of Nz ) and the maximum was set to be δmin z wavelength (a function of Nλ ). The spectral pitch was set to δmin λ (a function of Nλ ) in Eq. (16). In Fig. 4, each of the numerically calculated relative sparsities q, which assume a PSNR of 30 dB as the acceptable reconstruction fidelity, are shown with the blue dots. They are averages obtained using ten combinations of different mask patterns of the CA and distributions of the point sources. The 1 norm was chosen as the regularizer R in Eq. (6) for TwIST. The analytically calculated relative sparsities q shown with the red grid in Fig. 4 were from Eq. (13) with σ = 0.6, which was chosen experimentally. The numerical relative sparsity q approaches the analytical one, which is approximately independent of the numbers Nz and Nλ of depth channels and spectral channels in this case, as N x and Nλ increase, although there is a small difference between them based on the worst-case analysis. 6.
Conclusion
In summary, we extended SPICA to single-shot acquisition of multidimensional complex objects. The generalized system model and its reconstruction process were presented. As an example, assuming a depth-variant multispectral object, the system model and the theoretical conditions were derived and numerically demonstrated. In the multidimensional SPICA method, both the depth and spectrum dimensions are encoded by free space propagation, and the modulated depth and spectral dimensions are coherently and incoherently multiplexed, respectively. In the numerical demonstrations, we showed the effect of the error between the CA-based support constraint and the actual mask pattern in the reconstruction process. Other modeling errors, such as discretization [23], should also be investigated. The proposed imaging modality can be readily applied to other optical dimensions, such as time and polarization. Our approach consists of only a CA, a monochrome image sensor, and a coherent light source, and does not require any reference light. Thus, a compact and highly functional imaging system can be realized. In addition, the technique is applicable to incoherent holography and sensing in various spectral regions, such as X-ray and electron imaging.
#232104 - $15.00 USD © 2015 OSA
Received 8 Jan 2015; revised 13 Mar 2015; accepted 13 Mar 2015; published 7 Apr 2015 20 Apr 2015 | Vol. 23, No. 8 | DOI:10.1364/OE.23.009696 | OPTICS EXPRESS 9704
