Recent works have demonstrated that three-dimensional (3D) object reconstruction is possible from integral images captured in severely photon starved conditions. In this paper we propose an iterative approach to implement a maximum likelihood expectation maximization estimator with several types of regularization for 3D reconstruction from photon counting integral images. We show that the proposed algorithms outperform the previously reported approaches for photon counting 3D integral imaging reconstruction. To the best of our knowledge, this is the first report on using iterative statistical reconstruction techniques for 3D photon counting integral imaging.
©2011 Optical Society of America
Integral imaging [1–16] (InI) is an autostereoscopic three-dimension (3D) passive imaging technique originally proposed by G. Lippmann  a century ago and was revived in the last decade due to advances in digital imaging technology and optical and optoelectronic devices. Sensing of the 3D scene is carried out by capturing elemental images optically using a pickup microlens array and a detector array (Fig. 1 ) or by using multiple cameras. Integral imaging acquisition may be viewed as a multichannel system  where each channel generates an elemental image (EI). As such, there is a large redundancy in the captured data permitting visualization of occluded objects , 3D target recognition , and 3D image restoration in photon starved conditions [12–16]. 3D object reconstruction can be performed either optically or computationally. Optical reconstruction is depicted in Fig. 1(b). Computational reconstruction typically mimic the optical reconstruction, however has the additional flexibility of digitally manipulating the data to extract better visual information . In this paper we address the problem of computationally visualizing 3D objects from photon counting integral images. Following the work by Yeom et al.  demonstrating the possibility to recognize 3D object from InI with very small photon counts per elemental image, several approaches for extracting the visual information from InI captured in photon starved environments have been introduced [13–16]. In Ref.  Maximum Likelihood Estimation (MLE) with Poisson distribution was used for 3D passive sensing and visualization using InI for 3D image reconstruction· The photon counting detection model  was used to reconstruct 3D objects from photon counting images. Although MLE produced reasonably good results, simple MLE cannot provide optimal reconstruction because no prior information about the object is exploited. A follow up study  used Bayesian estimation with Gamma distribution as prior. This estimator is also not optimal since Gamma prior may not necessarily be the best way to describe object characteristics. Another study  considered truncated Poisson distribution which describes more accurately the photon counting conditions at extremely low photon, however it lacks an efficient way to include prior constrains.
In this work, we propose to reconstruct photon counting InI by using several penalized maximum likelihood expectation maximization (PMLEM) algorithms. In general, PMLEM algorithms effectively combines the Poisson statistics model associated with photon counting process, with penalties reflecting prior information or assumptions about the image (e.g., see  and references therein). MLEM algorithms were found to be very useful in various photon starved imagery, such as astronomy , nuclear medicine (e.g., Refs. [19–22]), and microscopy . Here we demonstrate its effectiveness in photon counting integral-imaging.
The paper is organized as following. In section 2 we briefly overview the 3D integral imaging; we present the PMLEM algorithm and explain its implementation in photon counting integral imaging. Then, in section 3, experimental results will be presented followed by concluding remarks in section 4.
2. Integral imaging with penalized maximum likelihood expectation maximization reconstruction
2.1 Reconstruction algorithms
Integral imaging (InI) is a three-dimensional (3D) imaging technique that can sense a 3D scene and reconstruct it as a 3D volumetric image. The InI acquisition scheme is depicted in Fig. 1. The reconstruction of the 3D object can be carried out optically or computationally [1,13]. Typically, computational reconstruction of volume pixels of the scene is accomplished by numerically simulating the optical reconstruction process; that is, by back propagating the captured light rays in Fig. 1(b). Image depth planes at arbitrary distances from the display microlens array are reconstructed by back propagating the elemental images (multi-view images) through a virtual pinhole array following ray optics [1,13]:
where is the reconstruction image at a specific depth plane and g represents the elemental images, and are the shifting factor between each lens at x and y directions respectively, K and L represent the number of elemental images at the CCD plane, k and l are indexes that pointing the pixel at such elemental image. For object located at distance z = from the lens of the imaging device, is the magnification factor.
Algebraically, relation (1) can be expressed as
where is a lexicographic vector representing backprojected reconstruction,, and g is a lexicographic vector representing the elemental images at the CCD plane. is the adjoin of the forward projection operator A defining the relation between the captured image g and the object f:Fig. 1). We note that back propagation model expressed by Eq. (1) holds in the geometrical approximation. In case that diffraction and aberrations are significant it can be easily generalized by accounting for local point spread function in A.
Equation (2) holds for the noiseless case. However, in low photon acquisition mode, the measurement vector g is subject to variations due to Poisson noise from the discrete nature of the arrival of the photon detection [13–16]. Since variance of a Poisson random variable equals to its mean, at low photon counts the noise added to g in Eq. (2) is significant. Our goal is to estimate the original 3D image f from g corrupted by Poisson noise. A well-known iterative solution for this problem is the PMLEM algorithm [18–23]:
where is the adjoin operator of representing the backward projection, β is a constant number controlling the amount of regularization induced by operator P. is the sensitivity vector with elements given by where ai,j is the i,jth element of the matrix A; hence is the sum of the jth object point response in the InI plane accounting for the number of EIs contributing to reconstruction of the point. The division in Eq. (4) should be taken as Hadamarad division (component-wise division).
which, for a well-conditioned A matrix, reduces to the well-known Richardson-Lucy deconvolution algorithm .
It can be seen that PMLEM algorithm in Eq. (4) uses the forward projection operator as well as the backward projection operator to estimate iteratively the scene based on the measurements. They can be seen as a set of successive projections and backprojections where the factor is the ratio between the captured image- g, and the one that would be obtained from the estimated object- at the kth iteration step. We notice that any positive initial guess yields a non-negative guess.
Since the acquisition is ill posed the reconstruction process needs to be regularized. This is done by including the term in the dominator of Eq. (4). This term penalizes the (k + 1)th guess according to the penalty operator P controlled by the parameter b. Penalty P is designed to reflect some prior assumption about the object, such as local smoothness of f(x,y). Consequently it is calculated on a patch surrounding the estimated pixel. The operator P is applied on the previous guess f (k), therefore the approach is called one-step-late .
For regularization of the MLEM algorithm applied for the photon counting integral imaging, we choose to investigate penalties implementing three types of priors: quadratic prior, median root prior, and total-variation prior.
2.2 PMLEM with quadratic prior
Probably the most direct prior for imposing smoothness on the image is the quadratic prior for which the penalty is given by 
where Si,j is a set of pixel in the neighborhood of pixel , and is a weighting function depending on the distance from the pixel.
2.3 Median root prior
Median root prior assumes that the image is locally smooth such that it passes unaltered through the median filter. This prior is introduced in the PMLEM by the nonlinear operator :
where is a median filter applied on a window Si,j in the neighborhood of pixel .
2.4 Total variation average penalty
Minimum total-variation is known as an effective prior for edge preserving regularization. An expression for the total variation penalty for PMLEM is derived in :22]. The total-variation constrain has been used extensively in image processing and it was previously proposed to penalizing the MLE algorithm for nuclear medicine [19,22] and image reconstruction in microscopy .
3. Experimental results
We applied the PMLEM algorithm on real and simulated integral images with the three penalties mentioned in the previous section. For the quadratic and median penalty we have Si,j to be a 3x3 window surrounding the pixel .
The PMLEM is an iterative algorithm; therefore it requires a stopping condition. The iterations are stopped when the means square error (MSE) between the estimated image at iteration k, , and the previous one is smaller than a threshold ,
where N and M are the number of pixels in the x and y coordinates respectively.
Table 1 show representative comparison PSNR results for the three types of penalties and the un-penalized MLEM. The PSNR is calculated as
where MSE is the mean square error between the reconstructed image from the photon starved image , and the reconstructed image from full photon statistics (i.e., regular illumination conditions) :
where i,j denote the pixel's index in x and y direction, respectively, and ROI is the region of interest area (i.e., car without background in Fig. 2 ). Table 1 reveals that total variations PLMEM (TV-PMLEM) lead to a higher PSNR of about 21.66dB. Figure 2 shows representative results obtained from real images captured with the InI setup described in Ref. . To evaluate the proposed 3D recognition system, we processed the same data presented in Ref. . obtained with the experiment described in the following. Elemental images of a toy car [Fig. 2(a)] were captured by moving a CCD camera laterally in both the x and y directions. The image sensor array had pixels and pixel size was. A white-light source was used to illuminate the object. By simulating the Poisson process  associated with a photon counting detector (e.g., single photon avalanche diodes) an array of elemental images was generated. The car was reconstructed at a distance of Z 0 = 370 mm from the sensor (see Fig. 2) using MLE and PMLEM algorithms. The initial guess f(0) in Eq. (4) was a positive-valued randomly generated image and parameter β in the range of 0.2-0.4 gave best results.
In Fig. 2 we present a representative comparison between the unpenalized MLEM to the PMLEM algorithms and the MLE algorithm. The reconstruction is carried out from only 100 photon counts per elemental image of 50635 pixels (i.e., average photon count of per pixel). Figure 2(a) shows the original toy car under usual illumination condition. Figure 2(b) shows an elemental image with approximately 100 counts. Figure 2(c) shows the reconstructed image using MLEM without penalty (Eq. (5)). In Fig. 2(d) we show the PMLEM reconstruction using quadratic prior. Figure 2(e) shows reconstruction with median prior penalty. It can be seen that, in concordance with PSNR results in Table 1, the median prior yield better visual reconstruction than the previous two. Best reconstruction results were obtained with TV-PMLEM, as demonstrated by Fig. 2(f). It can be seen that TV-PMLEM yields better reconstructions than other algorithms, both visually and in terms of PSNR. For comparison to previously published techniques, MLE reconstruction  is shown in Fig. 2(g). It can be seen that all the PMLEM algorithms give better results than the MLE algorithm. The bets improvement compared to MLE algorithm was obtained with TV-PMLEM algorithm; 2.98dB higher PSNR. To confirm the generality of this result we conducted a Monte Carlo simulative experiment involving 10 different objects. The simulations showed PSNR improvement of TV-PMLEM compared to MLE by 2.5 ± 0.6dB.
As demonstrated by Table 1 and Fig. 2, TV-PMLEM outperforms the other algorithms for very low photon counts that was considered in . Figure 3 shows the reconstruction PSNR for even lower photon counts for the InI in Fig. 2. It can be seen from Fig. 3 that TV-PMLEM produces similar results to MLE (PSNR≈17.7dB) with , that is, with only 1/4 of the number of photon counts per pixel.
In conclusion, we have shown that PMLEM outperforms some of the previously proposed methods such as (MLE) in photon counting integral imaging. Among the penalties considered total variation gave best results. TV-PMLEM produces same reconstruction results as previously proposed MLE with approximately 4 times less photons per pixel.
B. Javidi wishes to acknowledge support in part by DARPA and AFRL under contract number FA8650-09-C-7934. The opinions presented in this paper do not and are not mean to represent the opinions of DARPA and AFRL.
References and links
1. A. Stern and B. Javidi, “Three dimensional image sensing, visualization, and processing using integral imaging,” Proc. IEEE 94(3), 591–607 (2006). [CrossRef]
2. G. Lippmann, “La photographie integrale,” C. R. Acad. Sci. 146, 446–451 (1908).
4. B. Javidi, F. Okano, and J.-Y. Son, eds., Three-Dimensional Imaging, Visualization, and Display (Springer, 2008).
5. T. Okoshi, “Three-dimensional displays,” Proc. IEEE 68(5), 548–564 (1980). [CrossRef]
6. H. Hoshino, F. Okano, H. Isono, and I. Yuyama, “Analysis of resolution limitation of integral photography,” J. Opt. Soc. Am. A 15(8), 2059–2065 (1998). [CrossRef]
7. R. Martinez-Cuenca, G. Saavedra, M. Martinez-Corral, and B. Javidi, “Progress in 3-D multiperspective display by integral imaging,” Proc. IEEE 97(6), 1067–1077 (2009). [CrossRef]
9. F. Okano, J. Arai, K. Mitani, and M. Okui, “Real-time integral imaging based on extremely high resolution video system,” Proc. IEEE 94(3), 490–501 (2006). [CrossRef]
15. I. Moon and B. Javidi, “Three dimensional imaging and recognition using truncated photon counting model and parametric maximum likelihood estimator,” Opt. Express 17(18), 15709–15715 (2009). [CrossRef] [PubMed]
16. I. Moon and B. Javidi, “Three-dimensional recognition of photon-starved events using computational integral imaging and statistical sampling,” Opt. Lett. 34(6), 731–733 (2009). [CrossRef] [PubMed]
17. J. W. Goodman, Statistical Optics (Wiley, 1985).
18. M. Bertero, P. Boccacci, G. Desidera, and G. Vicidomini, “Image deblurring with Poisson data: from cells to galaxies,” Inverse Probl. 25(12), 123006 (2009). [CrossRef]
19. R. M. Lewitt and S. Matej, “Overview of methods for image reconstruction from projections in emission computed tomography,” Proc. IEEE 91(10), 1588–1611 (2003). [CrossRef]
20. S. Geman and D. Geman, “Stochastic relaxation, Gibbs distributions and Baysesian restoration of images,” IEEE Trans. Pattern Anal. Mach. Intell. PAMI-6(6), 721–741 (1984). [CrossRef]
22. V. Y. Panin, G. L. Zeng, and G. T. Gullberg, “Total variation regulated EM algorithm [SPECT reconstruction],” IEEE Trans. Nucl. Sci. 46(6), 2202–2210 (1999). [CrossRef]
23. N. Dey, L. Blanc-Feraud, C. Zimmer, P. Roux, Z. Kam, J. C. Olivo-Marin, and J. Zerubia, “Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution,” Microsc. Res. Tech. 69(4), 260–266 (2006). [CrossRef] [PubMed]