Open Access
Add to BibSonomyAbstract
We propose a preconditioning method to improve the convergence of iterative reconstruction algorithms in multiplexed imaging based on convolution-based compressive sensing with spatially coded point spread functions (PSFs). The system matrix is converted to improve the condition number with a preconditioner matrix. The preconditioner matrix is calculated by Tikhonov regularization in the frequency domain. The method was demonstrated with simulations and an experiment involving a range detection system with a grating based on the multiplexed imaging framework. The results of the demonstrations showed improved reconstruction fidelity by using the proposed preconditioning method.
©2011 Optical Society of America
1. Introduction
Single-shot acquisition of multi-dimensional or multi-channel optical data is an important hot topic in the optics field. Various imaging systems of this kind have been proposed [1, 2]. Those systems often used image sensors; however, image sensors can capture only a two-dimensional distribution of light intensity. Multi-shot imaging, multi-aperture imaging, or other smart optical coding schemes must to be employed to capture optical data of more than two dimensions with an image sensor. Compressive sensing (CS), which solves an inverse problem of an ill-posed linear system by employing a sparsity constraint [3–5], has been employed in single-shot acquisition of various multi-dimensional or multi-channel optical data, including range detection, spectral imaging, polarization imaging, and so on [6–10]. Also, CS with coded point spread function (PSF) convolution has been proposed [11–13] and applied to some multi-channel optical data acquisition techniques, including range detection, spectral imaging, and wide field-of-view imaging [6, 14, 15]. In the multi-channel data acquisition systems, each of the optical channels is convolved with spatially coded PSFs, and the encoded signals are multiplexed onto an image sensor. The object data can be reconstructed by CS algorithms. This framework can be applied to various kinds of optical information acquisition with a single shot [15].
To reconstruct the object data accurately, the correlations between the PSFs of the channels should be low, which means that the PSF must not be a delta function. This tends to increase the condition number of the system matrix, degrade the convergence property of iterative reconstruction algorithms, and increase the noise sensitivity. The problem is pronounced in incoherent optical systems because the PSF in the system must not be negative. Preconditioning methods have been proposed for some ill-conditioned and well-posed multi-shot imaging systems to boost the convergence in iterative reconstruction algorithms [16, 17]. The preconditioning methods employ Tikhonov regularization, which is useful for solving linear inversion problems, and the solutions have been used in various well-posed imaging systems [18, 19]. Here, to improve the convergence of the reconstruction, we modify the preconditioning methods in Refs. [16, 17] for the convolution-based CS framework for multi-channel data acquisition, which involves ill-posed linear inversion problems, and we demonstrate the effects.
2. The proposed preconditioning method
In the convolution-based CS framework for multi-channel data acquisition, optical signals on each channel in the object are convolved with an initially designed PSF, and the spatially encoded signals are multiplexed onto the detectors [15]. For simplicity, the system is represented by a one-dimensional model in this paper. The discretized system model can be written as
where 𝒢, 𝒫, and ℱ are the captured data, the spatial PSFs, and the object data, respectively, and x and c are indexes of the spatial coordinate and the channels.Equation (1) is rewritten in matrix notation as
where g ∈ ℝNx×1, Φ ∈ ℝNx×(Nx×Nc), and f ∈ ℝ(Nx×Nc)×1 are the vectorized captured data, the system matrix, and the vectorized object data, respectively; Nx and Nc are the numbers of detectors and channels of the object; ℝa × b is an a × b matrix of real numbers; C c ∈ ℝNx×Nx is a circulant matrix indicating the PSF for the c-th channel in the object; D c ∈ ℂNx×Nx, 0 ∈ ℝNx×Nx, and F ∈ ℂNx×Nx are a diagonal matrix indicating the optical transfer function (OTF) for the c-th channel in the object, a zero matrix, and a Fourier transform matrix, respectively; and ℂa×b is an a ×b matrix of complex numbers.The above system is ill-posed. To reconstruct the object, a CS method, which solves an ill-posed linear system by employing a sparsity constraint [3], is used. In this paper, a CS method solves the following problem:
where P, α, and ℛ(·) are a preconditioner matrix, a regularization parameter, and a regularizer, respectively; and ||·||ℓ 2 denotes ℓ 2 norm. Here, Tikhonov inversion of the system matrix Φ is chosen as the preconditioner matrix P to reduce the condition number of Φ [16]. In this case, PΦ is close to an identity matrix.The preconditioner matrix P ∈ ℝ(Nx×Nc)×Nx is given by
where is a diagonal matrix according to Eq. (4). The element in the i-th row and column in can be calculated with Tikhonov regularization by where d i ∈ ℂ1×Nc and are [(D 0)i (D 1)i ··· (D Nc−1)i] and . Here, (·)i denotes the element in the i-th row and column in a diagonal matrix; λ and I′ ∈ ℝNc×Nc are a regularization parameter and an identity matrix, respectively; and T and H show the transpose and Hermitian conjugate transpose matrices, respectively. A small λ realizes quick convergence, but on the other hand, it boosts noise on the captured data.3. Demonstrations
The proposed preconditioning method is demonstrated by simulations and an experiment. In the demonstrations, a CS method called the two-step iterative shrinkage/thresholding algorithm (TwIST) [20] is used to solve Eq. (5). TwIST is an iterative convex optimization algorithm using two previous estimates to improve the convergence properties. TwIST has often been used in recently demonstrated CS systems in the optics field. The number of iterations in TwIST was 2000. The regularization parameters of α in Eq. (5) and λ in Eq. (7) were chosen experimentally.
3.1. Simulations
Two three-dimensional objects f in the simulations are shown in Fig. 1. The size of the objects was 256 ×256 ×4. The object shown in Fig. 1(a) consisted of multiple rotated Shepp-Logan phantoms, which were sparse in the two-dimensional total variation (TV) domain [21]. The number of nonzero coefficients of the object in the two-dimensional TV domain was 8722. The other object shown in Fig. 1(b) was made from the House image in the Spectral image database [22]. The number of nonzero coefficients of the object in the three-dimensional discrete wavelet transform (DWT) domain was 5000.
Fig. 1 Objects in the simulations. (a) An object which is sparse in the two-dimensional TV domain and (b) an object which is sparse in the three-dimensional DWT domain.
Two sets of PSFs are shown in Fig. 2. The size of the PSF sets was 256 ×256 ×4. Each channel of the PSF set shown in Fig. 2(a) consisted of randomly located delta functions. The PSF set may be ideal to distinguish the channels’ data in the reconstruction, but the physical implementation is difficult. Each channel of the PSF set shown in Fig. 2(b) consisted of a single circular function. The diameters of the circular functions of the channels are different from each other. Depth from defocus or spectral imaging with chromatic aberration may be assumed in the PSF set. The captured data g of Fig. 1 with the PSFs in Fig. 2 is shown in Fig. 3. The size of the captured data was 256 ×256. The measurement signal-to-noise ratio (SNR) was 60 dB.
Fig. 2 Sets of PSFs. The sets are composed of (a) multiple delta functions and (b) single circular functions.
Fig. 3 The captured data of (a) Fig. 1(a) with the PSF in Fig. 2(a), (b) Fig. 1(a) with the PSF in Fig. 2(b), (c) Fig. 1(b) with the PSF in Fig. 2(a), and (d) Fig. 1(b) with the PSF in Fig. 2(b)
The reconstruction results with the original TwIST and the preconditioned TwIST of Fig. 1(a) are shown in Fig. 4. The peak SNRs (PSNRs) between Fig. 1(a) and Figs. 4(a), 4(b), 4(c), and 4(d) were 24.9 dB, 40.3 dB, 17.3 dB, and 25.8 dB, respectively. The reconstruction results with the original TwIST and the preconditioned TwIST of Fig. 1(b) are shown in Fig. 5. The PSNRs between Fig. 1(b) and Figs. 5(a), 5(b), 5(c), and 5(d) were 34.8 dB, 36.7 dB, 34.1 dB, and 34.9 dB, respectively. The reconstruction fidelity was improved by the proposed preconditioning method.
Fig. 4 The reconstructed data of Fig. 1(a) with (a) the original TwIST from Fig. 3(a), (b) the preconditioned TwIST from Fig. 3(a), (c) the original TwIST from Fig. 3(b), and (d) the preconditioned TwIST from Fig. 3(b).
Fig. 5 The reconstructed data of Fig. 1(b) with (a) the original TwIST from Fig. 3(c), (b) the preconditioned TwIST from Fig. 3(c), (c) the original TwIST from Fig. 3(d), and (d) the preconditioned TwIST from Fig. 3(d).
Figure 6 shows the relationship between the measurement SNR and the reconstruction PSNR obtained with the original TwIST and the preconditioned TwIST. The plots in Fig. 6(a), where the two-dimensional TV was used as the regularizer in Eq. (5), show that the preconditioning method is effective for high-SNR measurements. On the other hand, the plots in Fig. 6(b), where the three-dimensional DWT was used as the regularizer in Eq. (5), show that the preconditioning method is effective especially for low-SNR measurements. Figure 6 indicates improved reconstruction fidelity with a certain number of iterations by using the proposed preconditioning method. Note that the preconditioning method does not change the solution of the inversion problems. Therefore, the reconstruction results with infinite iterations in the original TwIST and the preconditioned TwIST are expected to be the same.
Fig. 6 Plots of reconstruction PSNRs from noisy measurements of (a) Fig. 1(a) and (b) Fig. 1(b) in the original TwIST and the preconditioned TwIST.
3.2. Experiment
The proposed preconditioning method was experimentally demonstrated with a range detection system based on the convolution-based CS framework for multi-channel data acquisition. The experimental setup is shown in Fig. 7. A phase grating (GTI50-03A manufactured by Thorlabs, line pitch: 300 lines/mm) was set in front of the imaging optics. The imaging optics was a COSMICAR television lens with a focal length of 16 mm. A monochrome detector array was used (Cool SNAP ES manufactured by Photometrics.) The number of pixels and the pixel pitch were 1392 ×1040 and 6.45 μm × 6.45 μm, respectively. The grating modulated the PSFs differently at each object distance, as shown in Fig. 7. Two planar objects, the printed characters “a” and “b”, were located at different distances from the imaging optics. The objects were passively illuminated with incoherent interior lights.
The captured data is shown in Fig. 8(a). The size was 193 ×446. The PSFs were two delta functions and the locations were varied with the object distance. The reconstruction results with the original TwIST and the preconditioned TwIST are shown in Figs. 8(b) and 8(c), respectively. The results show that the preconditioning method suppressed some artifacts in reconstruction.
Fig. 8 The experimental results. (a) the captured data and the reconstructed data with (b) the original TwIST and (c) the preconditioned TwIST.
4. Conclusions
We have proposed a preconditioning method to boost the convergence of iterative reconstruction algorithms in multiplexed imaging based on convolution-based CS with spatially coded PSFs. A preconditioner matrix using Tikhonov regularization in the frequency domain was introduced to reduce the condition number of the system matrix. The preconditioning method was demonstrated with simulations and an experiment. In the simulations, two PSF sets, one composed of randomly located delta functions and the other composed of single circular functions, were employed with TwIST using a sparsity constraint with a two-dimensional TV or a three-dimensional DWT. A range detection system with grating modulating PSFs at each distance was also experimentally demonstrated. The simulations and the experiment showed improved reconstruction fidelity with a certain number of iterations by using the preconditioning method.
References and links
1. A. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, 1988).
2. D. J. Brady, Optical imaging and spectroscopy (Wiley-OSA, 2009). [CrossRef]
3. D. L. Donoho, “Compressed sensing,” IEEE Trans. Info. Theory 52, 1289–1306 (2006). [CrossRef]
4. R. Baraniuk, “Compressive sensing,” IEEE Sig. Process. Mag. 24, 118–121 (2007). [CrossRef]
5. E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Sig. Process. Mag. 25, 21–30 (2008). [CrossRef]
6. D. J. Brady, K. Choi, D. L. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Express 17, 13040–13049 (2009). [CrossRef] [PubMed]
7. M. E. Gehm, R. John, D. J. Brady, R. M. Willett, and T. J. Schulz, “Single-shot compressive spectral imaging with a dual-disperser architecture,” Opt. Express 15, 14013–14027 (2007). [CrossRef] [PubMed]
8. R. Horisaki, K. Choi, J. Hahn, J. Tanida, and D. J. Brady, “Generalized sampling using a compound-eye imaging system for multi-dimensional object acquisition,” Opt. Express 18, 19367–19378 (2010). [CrossRef] [PubMed]
9. M. Shankar, N. P. Pitsianis, and D. J. Brady, “Compressive video sensors using multichannel imagers,” Appl. Opt. 49, B9–B17 (2010). [CrossRef] [PubMed]
10. A. Ashok and M. A. Neifeld, “Compressive light field imaging,” Proc. SPIE 7690, 76900Q (2010). [CrossRef]
11. J. Romberg, “Compressive sensing by random convolution,” SIAM J. Imaging Sci. 2, 1098–1128 (2009). [CrossRef]
12. R. F. Marcia and R. M. Willett, “Compressive coded aperture superresolution image reconstruction,” in “IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2008),” (2008), pp. 833–836.
13. Y. Rivenson, A. Stern, and B. Javidi, “Single exposure super-resolution compressive imaging by double phase encoding,” Opt. Express 18, 15094–15103 (2010). [CrossRef] [PubMed]
14. J. Hahn, S. Lim, K. Choi, R. Horisaki, and D. J. Brady, “Video-rate compressive holographic microscopic tomography,” Opt. Express 19, 7289–7298 (2011). [CrossRef] [PubMed]
15. R. Horisaki and J. Tanida, “Multi-channel data acquisition using multiplexed imaging with spatial encoding,” Opt. Express 18, 23041–23053 (2010). [CrossRef] [PubMed]
16. N. Nguyen, P. Milanfar, S. Member, and G. Golub, “A computationally efficient superresolution image reconstruction algorithm,” IEEE Trans. Image Proc. 10, 573–583 (2001). [CrossRef]
17. K. Choi, R. Horisaki, J. Hahn, S. Lim, D. L. Marks, T. J. Schulz, and D. J. Brady, “Compressive holography of diffuse objects,” Appl. Opt. 49, H1–H10 (2010). [CrossRef] [PubMed]
18. A. Ashok and M. A. Neifeld, “Pseudorandom phase masks for superresolution imaging from subpixel shifting,” Appl. Opt. 46, 2256–2268 (2007). [CrossRef] [PubMed]
19. A. Mahalanobis, M. Neifeld, V. K. Bhagavatula, T. Haberfelde, and D. Brady, “Off-axis sparse aperture imaging using phase optimization techniques for application in wide-area imaging systems,” Appl. Opt. 48, 5212–5224 (2009). [CrossRef] [PubMed]
20. 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). [CrossRef]
21. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Phys. D 60, 259–268 (1992). [CrossRef]
22. “Spectral image database,” http://spectral.joensuu.fi/multispectral/spectralimages.php.
