1. Introduction

Integral imaging [1–16] (InI) is an autostereoscopic three-dimension (3D) passive imaging technique originally proposed by G. Lippmann [2] 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 [1] where each channel generates an elemental image (EI). As such, there is a large redundancy in the captured data permitting visualization of occluded objects [3], 3D target recognition [1], 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 [1]. In this paper we address the problem of computationally visualizing 3D objects from photon counting integral images. Following the work by Yeom et al. [12] 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. [13] 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 [17] 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 [14] 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 [15] 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.

 

Fig. 1 (a) forward projection of a 3D object in the InI plane through microlens array. (b) Backward projection.

Download Full Size | PDF

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 [18] and references therein). MLEM algorithms were found to be very useful in various photon starved imagery, such as astronomy [18], nuclear medicine (e.g., Refs. [19–22]), and microscopy [23]. 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 $f^{BP}$ is the reconstruction image at a specific depth plane and g represents the elemental images, $S_x$ and $S_y$ 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 = $z_0$from the lens of the imaging device, $M_0$ is the magnification factor.

Algebraically, relation (1) can be expressed as

where $f^{BP}$is a lexicographic vector representing backprojected re­con­struc­tion,$f^{BP}=vect\left(f\left[n,m\right]\right)$, and g is a lexicographic vector representing the elemental images at the CCD plane. $A^T$ is the adjoin of the forward projection operator A defining the relation between the captured image g and the object f:

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 $A^T$ is the adjoin operator of representing the backward projection, β is a constant number controlling the amount of regularization induced by operator P. $s_j$ is the sensitivity vector with elements given by $s_j={\displaystyle \sum _{i=1}^m{a}_{i,j}},j=1,\mathrm{...}n$ where a_i,j is the i,jth element of the matrix A; hence $s_j$ 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).

It is noted that in the case β = 0 Eq. (4) represents the MLEM [18] algorithm:

which, for a well-conditioned A matrix, reduces to the well-known Richardson-Lucy deconvolution algorithm [18].

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$\frac{g}{A\cdot {\widehat{f}}^{(k)}}$ is the ratio between the captured image- g, and the one that would be obtained from the estimated object- ${\widehat{g}}^{(k)}=A\cdot {\widehat{f}}^{(k)}$ at the kth iteration step. We notice that any positive initial guess$f^{(0)}$ yields a non-negative guess${\widehat{f}}^{(k)}$.

Since the acquisition is ill posed the reconstruction process needs to be regularized. This is done by including the term $\beta \cdot P{\widehat{f}}^{\left(k\right)}$ 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$f_{i,j}^{(k)}$. The operator P is applied on the previous guess f ^(k), therefore the approach is called one-step-late [19].

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 [20]

where S_i,j is a set of pixel in the neighborhood of pixel ${\widehat{f}}^{\left(k\right)}[i,j]$, and $w_{i,j}$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 [21]:

where $Med\left\{f[i,j];S_{i,j}\right\}$ is a median filter applied on a window S_i,j in the neighborhood of pixel $f[i,j]$.

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]:

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 S_i,j to be a 3x3 window surrounding the pixel $f[i,j]$.

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, $f^{(k)}$, and the previous one $f^{(k-1)}$ is smaller than a threshold ${\epsilon }_{th}$,

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

Table 1. Reconstruction results with maximum likelihood estimation (MLE), maximum likelihood expectation maximization (MLEM) and penalized MLEM (PMLEM) with three types of penalties

View Table

where MSE is the mean square error between the reconstructed image from the photon starved image ${\widehat{f}}_{}^{est}$, and the reconstructed image from full photon statistics (i.e., regular illumination conditions) ${\widehat{f}}_{}^{full}$:

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. [16]. To evaluate the proposed 3D recognition system, we processed the same data presented in Ref. [16]. 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 $2028\times 2044$ pixels and pixel size was$9\mu m\times 9\mu m$. A white-light source was used to illuminate the object. By simulating the Poisson process [15] associated with a photon counting detector (e.g., single photon avalanche diodes) an array of $10\times 10$ elemental images was generated. The car was reconstructed at a distance of Z ₀ = 370 mm from the sensor (see Fig. 2) using MLE and PMLEM algorithms. The initial guess f⁽⁰⁾ in Eq. (4) was a positive-valued randomly generated image and parameter β in the range of 0.2-0.4 gave best results.

 

Fig. 2 Reconstruction at Z₀ = 370mm using InI with $\widetilde{Np}\approx 2\times {10}^{-3}$ photon/pixel, using: (a) normal illuminations conditions, (b) one elemental image, (c) MLEM without penalty (5), (d) PMLEM with prior reflecting quadratic penalty, (e) PMLEM with median penalty, (f) PMLEM with total-variation (g) using MLE.

Download Full Size | PDF

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 $\widetilde{Np}\approx 2\times {10}^{-3}$ 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 [13] 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 [16]. 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 $\widetilde{Np}\approx 0.5\times {10}^{-3}$, that is, with only 1/4 of the number of photon counts per pixel.

 

Fig. 3 Reconstruction PSNR for very low photon counts.

Download Full Size | PDF

4. Conclusions

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.