1. Introduction

A compound-eye imaging system is a promising computational imaging modality. Compound-eye optics have enabled light-field acquisition[1] and device compactness[2, 3]. Thin observation module by bound optics (TOMBO) is a representative example of a compound-eye imaging system[4].

An advantage of compound-eye imaging systems is that they permit diverse data acquisition schemes. Different lenslets may create different encodings. For example, time detection based on the encoding concept has been proposed[5], and range detection in [6] can be considered as a system based on the concept. These compact systems reconstruct a three-dimensional object from a two-dimensional measurement where the size is the same as that of an axial plane of the object.

This paper proposes generalized sampling approaches for multi-dimensional object acquisition using TOMBO. In the proposed system, an object is acquired with coding and multiplexing in a two-dimensional snapshot. In particular, the coding schemes in [5, 6] are extended for multi-dimensional data acquisition of various objects.

There could be multiple choices for coding schemes for multi-dimensional object acquisition such as coded aperture imaging and multi-shot imaging. These schemes differ by design constraints. This paper considers the compactness of hardware and single-shot object acquisition capability as critical design constraints. As indicated by a large literature, TOMBO imaging modality is one of such techniques that can implement a compact system that can meet our design constraints. This motivates us to investigate its potential as a compressive imaging technique.

In this paper, the mathematical model of the proposed system and examples of the coding schemes for spectral and polarization imaging techniques are presented. Simulation results of the proposed system are shown. The implementations are inspired by [7, 8, 9]. The previously presented systems have a tradeoff between the spatial and axial resolutions. For example, in [7, 8], the number of the spectral or polarization channels is roughly proportional to that of the lenses. Increasing the number of lenses reduces the spatial resolution. The approaches proposed in this paper may compensate for the tradeoff by leveraging compressive sampling[10].

A constrained optimization technique to incorporate sparsity in some basis of an object estimate is used for reconstruction. The reconstruction method is inspired by compressive sampling [10]. In compressive sampling, the systems should satisfy some assumptions stated in section2 for accurate reconstruction. The proposed system is compared to a theoretical baseline sensing system which is a Gaussian random sensing matrix. Several systems based on sparse reconstruction have been demonstrated and have shown promising results[11, 12, 13].

 

Fig. 1. Cross section view of TOMBO. ν, O_u, and L_u are the spatial dimension, the center position, and the position of a lenslet in the u-th unit, respectively.

Download Full Size | PDF

Section 2 provides a brief background on TOMBO and compressive sampling. Section 3 describes a general model for multi-dimensional TOMBO imaging. Section 4 presents examples of coding schemes. Simulation results are given in section 5.

2. Background

2.1. TOMBO

In a simplified conceptual model, TOMBO consists of lenslets and a detector array as shown in Fig. 1. An imaging structure associated with a lenslet is called a unit[4]. Each unit produces a low-resolution (LR) image on the detector array.

When the number of units is N_u × N_u in a square arrangement, the focal-length and the diameter of the lenslet need to be N_u times smaller than those of the corresponding conventional full-aperture system to obtain the same field of views. This results in a LR image whose size is N_u times smaller than that of an image produced by the full-aperture system. The thickness and depth-of-field of TOMBO are N_u times shorter and N_u ² times longer, respectively. This allows for compact hardware with a large depth-of-field. Objects are often assumed to be located within the depth-of-field, and the lenslets are assumed to be aberration-free[4, 14]. These assumptions are made throughout this paper, unless otherwise stated.

2.2. Compressive sampling

The proposed system model in this paper forms an underdetermined linear system of equations as described in section 3. Compressive sampling (CS) is a theoretical framework for solving an underdetermined system[10, 15]. The reconstruction method in this paper is inspired by CS.

A linear system model can be written as

where $g\in {ℝ}^{N_g\times 1}$ , $\mathbf{\Phi }\in {ℝ}^{N_g\times N_f}$ , $f\in {ℝ}^{N_f\times 1}$ , $\mathbf{\Psi }\in {ℝ}^{N_f\times N_{f\prime }}$ , and $\beta \in {ℝ}^{N_f^{\prime }\times 1}$ are a measurement vector, a sensing matrix, an object vector, a basis matrix, and a transform coefficient vector, respectively. ${ℝ}^{N_i\times N_j}$ denotes a N_i × N_j matrix of real numbers. We consider the case where N_g ≪ N_f.

Let s denote the number of non-zero coefficients in β. CS indicates that, for accurate reconstruction, Θ should satisfy a sufficient condition for any s-sparse β. The condition is called the restricted isometric property (RIP) defined by,

where c_s ∈ (0,1) is a constant and ||·||² ₂ denotes an ℓ₂-norm[16]. Λ is a subset of indices supporting s nonzero coefficients in β. β _Λ and Θ _Λ are elements of β and columns of Θ that support the s coefficients. If c_s is close to 0, Eq. (2) indicates that Θ _Λ preserves the Euclidean length of ${\beta }_{\Lambda }.\mu (\mathbf{\Phi },\mathbf{\Psi })\in [1,\sqrt{N_f\prime }]$ defined as

is called coherence. Φ(i, :), Ψ(:, j), and 〈,〉 are the i-th row of Φ, the j-th column of Ψ, and an inner product, respectively. When the coherence is small, Φ and Ψ are said to be incoherent. The number of measurement components required for accurate reconstruction is given as

where c is a constant[15]. According to CS theory[15], if Θ satisfies RIP (Eq. (2)), measurements are with high probability sufficient to accurately estimate β. An accurate estimate of the s nonzero coefficients in β can be obtained by solving

where || · ||₁ denotes ℓ₁ norm.

3. A mathmatical model for proposed acquisition schemes

Let F (x,y,z ₀, ⋯, $z_{N_n-1}$ denote a continuous density function representing a multi-dimensional object. x and y represent spatial dimensions, and z ₀, ⋯, $z_{N_n-1}$ represent the other dimensions dependent on the application. x = 0 and y = 0 are defined as the center of a detector array. For simplicity, the y dimension is omitted. Extending to higher dimensions may be readily achieved with small modifications of the model.

3.1. Continuous model

In the proposed system, a multi-dimensional object is integrated on to detectors with one of two coding schemes, as demonstrated in Fig. 2. In one of the coded integrations inspired by [6], an object is sheared by an optical element, and the sheared optical signal is integrated on to a detector array. In the shear-transformation, each axial plane in an object is shifted along the x axis as shown in Fig. 2(a). In [6], the shift corresponds to a parallax. In another coded integration inspired by [5], an object is multiplied with a weight distribution, and the weighted optical signal is integrated on to a detector array as shown in Fig. 2(b). The weight distribution is a continuous function of z. In [5], the weight distribution corresponds to an exposure time. The two schemes are referred to as sheared integration (SI) and weighted integration (WI), respectively.

We denote integrated data associated with the u-th unit as G_u(ν), where ν denotes the spatial dimension in a unit as shown in Fig. 1. ν in the u-th unit is defined as ν = x − O_u, where O_u is the center of the u-th unit. G_u(ν) is expressed as

where S_n,u(z_n) and W_n,u(z_n) show a shift in SI and a weight distribution in WI of the z_n dimension in the u-th unit, respectively. L_u is the center position of the u-th lenslet on the ν axis as shown in Fig. 1. For simplicity, N_n = 1 is assumed, and a subscript n is omitted. Eq. (6) can be rewritten as

 

Fig. 2. Coding schemes in TOMBO. (a) Sheared integration and (b) weighted integration in a unit.

Download Full Size | PDF

3.2. Discretization model

A discrete object $\in {ℝ}^{N_x\times N_z}$ and a discrete integrated data $\in {ℝ}^{N_x\times N_u}$ are denoted by F̃ (l,m) = F(lΔ_x,mΔ_z) and G̃′_u(i) = G_u(iΔ_x) using notation similar to that in [17], where a tilde indicates a discrete data. G̃′ is an intermediate data before sampling by the detectors. l, m, and i are integer variables of the x, z, and ν axes in a discretization model, respectively. Δ_x and Δ_z are the pixel pitches along the x and z axes in a discrete object. G̃′_u(i) is sampled by detectors. The measurement data $\in {ℝ}^{N_v\times N_u}$ is expressed as G̃_u(j) = ∑_i G̃′_u(i) rect((iΔ_x − jΔ_ν − D_u)/Δ_ν), where N_ν, j, Δ_ν, and D_u are the number of detectors in a unit, an index for the detectors in a unit, the pixel pitch of the detectors, and the center of the center detector on the ν axis in the u-th unit. Then, the measurement data G̃ can be written as

where S̄_u(m) = ⌊(L_u + S_u(mΔ_z) − D_u)/Δ_x + 0.5⌊ and W̄_u(m) = W_u(mΔ_z)Δ_z· ⌊·⌋ is the floor function.

3.3. System matrix

We assume that Δ_x = Δ_ν/N_u and N_x = N_ν N_u, which are both typical assumptions in TOMBO imaging [4, 5, 6, 14]. Thus, the numbers of elements in the measurement data and the object are N_g = N_x and N_f = N_x N_z.

As indicated in Eq. (8), the m-th axial plane in F̃ is shifted by S̄_u(m) and multiplied with W̄_u(m). $C{\prime }_{m,u}\in {ℝ}^{N_x\times N_x}$ which denotes the coding operation for the m-th axial plane in F̃ in the u-th unit is expressed as,

where C′ _m,u (p,q) is the (p,q)-th element in the matrix C′ _m,u.

$C_u\in {ℝ}^{N_f\times N_f}$ represents the coding operation implemented by the TOMBO system on the object f (F̃ in vector form) in the u-th unit and is written as

where O is a N_x × N_x zero matrix.

The matrix $Q\in {ℝ}^{N_x\times N_f}$ which sums all of the axial layers is defined by

where $I\in {ℝ}^{N_x\times N_x}$ denotes an identity matrix.

The downsampling matrix $T\in {ℝ}^{N_v\times N_x}$ can be defined by

where 1 and 0 denote a N_u × 1 vector whose elements are all 1 and a N_u × 1 vector whose elements are all 0, respectively. A superscript T indicates a transpose of a matrix.

Therefore, the measurement data G̃ on the u-th unit is TQC_uf, which means, firstly, an object is coded in each unit by C_u, secondly, the coded data is integrated on to a detector array by Q, and lastly, the integrated data is downsampled with detectors by T. The sensing matrix $\mathbf{\Phi }\in {ℝ}^{N_g\times N_f}$ is expressed by

4. Implementation of proposed acquisition schemes

The proposed coding schemes can be implemented for a wide array of practical applications. Each application would rely on some physical optical elements to implement the coding scheme expressed by Eqs. (6) or (7). In this section, we present examples of the coding scheme for spectral imaging and polarization imaging. Using similar schemes, physical coding strategies for range, time, spectrum, polarization, large dynamic range, and wide field-of-view may be available.

Physical codings for spectral imaging using SI and WI are illustrated in Fig. 3. SI for spectral imaging can be implemented by using dispersive elements (e.g., prisms). The elements in each unit have different dispersion directions as shown in Fig. 3(a). The dispersion results in different shifts for each spectral slice. In Eq. (7), z represents the wavelength. The shift corresponds to S_u(z).

WI for spectral imaging may be implemented with multi-band pass filters placed above or below the lenslet as shown in Fig. 3(b). Each of the filters has different pass-bands. Pass-bands and stop-bands are represented with W_u(z) = 1 and W_u(z) = 0 in Eq. (7), respectively. A stack of bandstop filters or a patch of bandpass filters may be used to substitute for the multi-band pass filter.

 

Fig. 3. Cross section views of TOMBO for spectral imaging with (a) SI and (b) WI.

Download Full Size | PDF

 

Fig. 4. Top views of TOMBO for polarization imaging with (a) SI and (b) WI. Arrows, dots, circles, and shaded areas indicate directions of polarization, centers of shifted images, lenslets, and polarization plates, respectively.

Download Full Size | PDF

Figure 4 shows a conceptual diagram for polarization imaging with the proposed codings. SI for polarization imaging may be performed with birefringent linear polarizers[18]. The elements split an incident ray into two polarized rays. Hence, an image at each polarization angle is shifted. Each unit has different shift for each polarization angle as shown in Fig. 4(a). Here, z represents a linear polarization angle. The shift corresponds to S_u(z) in Eq. (7).

WI for polarization imaging may be performed with polarization plates. Polarization plates with different linear polarization angle are placed above or below the lenslet as shown in Fig. 4(b). The weight distribution is expressed as W_u(z) = cos²(P_u − z)[18], where P_u is the polarization angle in the u-th unit. A patch of polarization plates, where each plate has a different polarization angle, allows flexibility in the design of a weight distribution.

5. Simulation of the proposed concept

The concept of multi-dimensional TOMBO imaging was verified through application independent simulations. These general simulations could readily modified for a specific application like those mentioned in the previous section.

A method called two-step iterative shrinkage/thresholding algorithm (TwIST) [19] was used for reconstruction. TwIST is an interative convex optimization algorithm that uses two previous estimates to improve convergence properties for the problem described by Eq. (5).

For simplicity, a shift in SI was assumed as S_u(z) = (A_uz + B_u)Δ_x. A_u and B_u are a gradient and a bias, respectively, of the shear-transformation in the u-th unit defined as A_u = (−2u/(N_u − 1) +1)A ₀ and B_u = −A_uN_zΔ_z/2. For example, A ₀ = 1.0 of N_u = 3 indicates that A ₀ = 1.0, A ₁ = 0.0, and A ₂ = −1.0. A shift at the center axial plane, where z = N_zΔ_z/2, is S_u(z) = 0.0. A_u and B_u of the y axis is the same as those of the x axis. A weight distribution in WI was assumed to be a binary pattern. In the m-th axial plane, h units were set as W_u(mΔ_z) = 1 in Eq. (7). The h units were randomly chosen, while the other N_u ² − h units were set as W_u(mΔ_z) = 0. In this case, the maximum number of separable axial planes is fixed ${}_{{}_{N_u}2}{C}_h$ . A lenslet’s position L_u in Eq. (8) was randomly set in each unit. The range was [−Δ_ν/2, Δ_ν/2], where Δ_ν is the pixel pitch of the detectors. The position of the center detector in a unit is D_u = 0.

Figure 5 shows a simulation of four-dimensional data acquisition using the two TOMBO coding schemes. An object whose size is 128 × 128 × 4 × 2 and measurement data whose size is 128 × 128 are shown in Figs. 5(a) and 5(b). The compression ratio is 8, which is calculated as N_f/N_g, where N_f and N_g are the numbers of elements in an object and measurement data, respectively. In Fig. 5, the object and the simulation results are reshaped to 128 × 128 × 8 for display. SI with A ₀ = 1.0 and WI with h = 3 were used for the z ₀ and z ₁ axes, respectively. The measurement signal-to-noise ratio (SNR) in the presence of additive white Gaussian noise and the number of units were 30 dB and 2 × 2, respectively. The object estimate sparsity in gradients was enforced using the total variation (TV) [20]. Two-dimensional TV was applied independently for each axial plane as ${\Sigma }_{l_x}{\Sigma }_{l_y}{\Sigma }_{m_0}{\Sigma }_{m_1}\mid \nabla {\left[\tilde{F}(l_x,l_y,m_0,m_1)\right]}_{l_x,l_y}\mid$ , where l_x and l_y are indices for the x and y axes in a discrete object. ∇[·]l_x,l_y is a two-dimensional gradient vector for the x and y directions, and ∣ · ∣ denotes the magnitude of the gradient vector. The object consists of multiple Shepp-Logan phantoms, which is sparse in two-dimensional TV domain. The total number of non-zero gradient values was s = 3242. The reconstruction results with TwIST and the Richardson Lucy method (RL)[21, 22] are compared in Figs. 5(c) and 5(d). Their peak signal-to-noise ratios (PSNR) were 32.1 dB and 19.4 dB, respectively. The PSNR is found by computing 20log₁₀( $\frac{\mathrm{MAX}}{\sqrt{\mathrm{MSE}}}$ ), where MAX and MSE represent the maximum of the signal values and the mean squared error between two signals, respectively[23].

CS object reconstruction accuracy may be estimated using a correlation between columns of Θ, which is the multiplication of a sensing matrix Φ and a basis matrix Ψ, in Eq. (1)[24]. When a correlation between two columns in Θ is high, it is difficult to resolve the two components in a transform coefficient vector β corresponding to the columns in Θ. So that, the reconstruction accuracy depends on not only Φ but also Ψ.

When a two-dimensional basis is used for each axial plane as in the previous simulation, the reconstruction accuracy along the axial direction in an object estimate may be roughly predicted based on the correlation between columns of Φ corresponding to two axial planes. Let ${\varphi }_m\in {ℝ}^{N_x\times N_x}$ denote

From Eqs. (1), (13), (14), and the assumption to use a two-dimensional basis, a sensing matrix, a basis matrix, and Θ can be rewritten as

respectively, where $\psi \in {ℝ}^{N_x\times N_x}$ and $O\in {ℝ}^{N_x\times N_x}$ are a two-dimensional basis matrix for each axial plane and a N_x × N_x zero matrix. If the correlation between a column in ϕ_m and that on another axial plane is high, then the corresponding correlation between a column in ϕ_mψ and that on another axial plane may be high. In this case, it is difficult to resolve the axial planes. When ∣A ₀∣ is small or h is large, the correlation between a column in ϕ_m and that on another axial plane is high. For example, Fig. 5(e) shows a reconstruction result where SI with A ₀ = 0.2 was used for the z ₀ axis. The reconstruction accuracy along the axial direction with A ₀ = 0.2 was lower than that with A ₀ = 1.0.

Figure 6 shows another simulation of five-dimensional data acquisition with the discrete wavelet transform (DWT). The sizes of the object in Fig. 6(a) and the measurement data in Fig. 6(b) were 128 × 128 × 2 × 2 × 2 and 128 × 128. The compression ratio is 8. Two-dimensional DWT was applied for each axial plane. The object consists of multiple natural images where the small coefficients in two-dimensional DWT were truncated. In two-dimensional DWT domain, the total number of non-zero DWT coefficients across all the planes was s = 2000. SI with A ₀ = 3.0, WI with h = 12, and WI with h = 12 were used for the z ₀, z ₁, and z ₂ axes, respectively. The measurement SNR and the number of the units were 30 dB and 4 × 4, respectively. The reconstruction results with TwIST and RL are compared in Figs. 6(c) and 6(d). Their _PSNRs were 24.5 dB and 15.4 dB, respectively. Figure 6(e) shows a reconstruction result where WI with h = 15 was used for the z ₁ and z ₂ axes. The reconstruction accuracy along the axial direction with h = 15 was lower than that with h = 12.

Figure 7 illustrates the sensitivity of the reconstructions to noise as represented by a curve relating the measurement SNR to the reconstruction PSNRs. Also, the performance is compared to that of an ideal Gaussian random compressive sensing matrix, which is known to require a (optimally) small number of measurements to satisfy the RIP compared to what the proposed systems would require. Since the proposed systems usually have a worse RIP meaning that more measurements are required to obtain higher reconstruction accuracy, they present worse reconstruction accuracy compared to that of the Gaussian random matrix. However, such random sensing matrices would be very difficult to physically implement in general with the current technology. In addition, it is not clear how such random sensing systems may provide the compactness of physical systems and snapshot acquisition functionality, which are benefits of the proposed approach.

6. Conclusions

We proposed a generalized sampling approach for multi-dimensional object acquisition using TOMBO. The sampling uses multi-dimensional sheared and weighted integration in each unit. The mathematical model and some examples of the proposed measurement approach were presented. The simulation demonstrated reconstruction of an object with the number of elements totaling eight times that of the measurement data. A method inspired by compressive sampling was used in the reconstruction. These schemes enable us to acquire a multi-dimensional object with a single two-dimensional measurement by a compact imaging system. Also, these schemes extend abilities of compound-eye imaging systems to various applications.

 

Fig. 5. Simulation results with total variation. (a) A four-dimensional object (∈ ℝ^{128×128×4×2}), where indices of axial planes are shown under each axial plane, (b) a measurement data, (c) a reconstruction with TwIST, (d) a reconstruction with RL, and (e) a reconstruction with TwIST using a small ∣A ₀∣.

Download Full Size | PDF

 

Fig. 6. Simulation results with discrete wavelet transform. (a) A five-dimensional object (∈ ℝ^{128×128×2×2×2}), (b) a measurement data, (c) a reconstruction with TwIST, (d) a reconstruction with RL, and (e) a reconstruction with TwIST using a large h.

Download Full Size | PDF

 

Fig. 7. Plots of reconstruction PSNR from noisy measurements in the proposed system and a baseline sensing system which is a Gaussian random sensing matrix. (a) Plots with the object, the parameters, and the basis used in Fig. 5(c) and (b) plots with the object, the parameters, and the basis used in Fig. 6(c).

Download Full Size | PDF

A useful avenue for future work is to analyze theoretical properties of the proposed systems. It would be interesting to see how many more measurements would be required in general for the proposed systems to produce a certain accuracy, which is related to the validity of the sparsity assumption in the proposed systems. Also, it would be very useful to find a more efficient sparsity transformation that provides a better RIP and a better sparse representation of the objects of interest. Furthermore, we plan to investigate other coding schemes that may provide a better RIP overall to better exploit the sparsity assumption.