Compressive spectral imaging system based on liquid crystal tunable filter

Xi Wang; Yuhan Zhang; Xu Ma; Tingfa Xu; Gonzalo R. Arce

doi:10.1364/OE.26.025226

1. Introduction

Spectral imaging systems based on liquid crystal tunable filters (LCTF) have been widely used in remote sensing, biomedicine, food industry and other fields due to their portability, fast tunability, convenient controllability, good image quality and low cost [1–7]. LCTF-based spectral imaging systems capture the three-dimensional data cube of the scene, including two spatial dimensions and one spectral dimension. As shown in Fig. 1(a), a typical LCTF-based hyperspectral imaging system consists of an imaging lens, followed by an LCTF and a detector array [2]. In the data acquisition stage, only a narrow-band spectrum of the scene can pass through the LCTF, and a quasi-monochromatic image of the scene is then acquired by the detector array at a time. In order to obtain the hyperspectral data cube across multitudes of wavelengths, the scene is spectrally scanned by tuning the center wavelength of the LCTF.

Fig. 1 (a) Traditional LCTF-based spectral imaging system, and (b) the sketch of the proposed CATF spectral imaging system.

Download Full Size | PDF

Traditionally, the LCTF was considered an ideal spectral filter with an approximately impulsive transmission function, and the filter output is approximated to a monochromatic image corresponding to its center wavelength [8]. However, practical LCTFs always have broader than desired bandwidths, and their transmittance is generally depicted using the real spectral signatures of LCTFs. Thus, the spectral resolution of the imaging system is limited by the bandwidth of the LCTF. Most of current LCTFs are designed based on Lyot filters [9]. This type of LCTFs can obtain narrower bandwidth by increasing the stage of LCTF or inserting additional bandpass filters [10–12]. However, the energy transmittance of light reduces with narrowing bandwidth, thus leading to inferior image quality and longer image acquisition time. Therefore, a shortcoming of LCTF-based spectral imagers is the inherent trade-off between the spectral resolution and optical throughput [2, 11]. On the other hand, the spatial resolution of traditional LCTF-based spectral imaging systems is essentially limited by the resolution of detector. However, high-resolution detector arrays are costly and difficult to manufacture, especially for the infrared waveband [13]. Another drawback of the traditional LCTF-based spectral imagers is that they need to densely scan and sample along the spectral dimension to obtain the hyperspectral images. The hyperspectral images include abundant and highly compressible data that is often compressed before storage and transmission [14–16].

In order to overcome these limitations, this paper develops a super-resolution LCTF-based hyperspectral imaging system using compressive sensing (CS) approaches. CS theory provides a way to acquire a high-dimensional signal with much lower sampling rate than that required by the Nyquist sampling theorem [17–19]. In the past, several compressive spectral imaging approaches have been proposed, such as the coded aperture snapshot spectral imager (CASSI) [20–23], the compressive sensing hyperspectral imager [24], the compressive hyperspectral imaging by spectral and spatial operators (CHISSS) [25] and so on [26]. However, the aforementioned compressive hyperspectral imaging systems are not LCTF-based, but use dispersive elements to discriminate different spectral components. Recently, August, et al. proposed an approach using the liquid crystal phase retarder to modulate and compress the spectra of hyperspectral images [27, 28]. However, this approach is inadequate to compress the hyperspectral data in spatial domain, and thus inherently limiting the compression ratio.

In contrast, the proposed spectral imager in this paper is capable of compressing the hyper-spectral data in both spatial and spectral domains simultaneously. We refer to the proposed system as the coded aperture tunable filter (CATF) spectral imager. The sketch of CATF spectral imager is illustrated in Fig. 1(b). The accurate transmission functions are used to represent the spectral response characteristics of the practical LCTFs. For a given center wavelength, the LCTF collects spectral components around the center wavelength, and modulates the spectrum by the amplitudes of the transmission functions. Then, the filtered spectral images are spatially modulated by a coded aperture with higher resolution than the detector. The coded aperture is implemented by the digital micromirror device (DMD), having an array of micromirrors that can be independently adjusted to create different block-unblock coded patterns. Finally, the spectral images are projected and multiplexed onto a low-resolution monochrome complementary metal oxide semiconductor (CMOS) detector via the relay lens. Each pixel on the detector will receive the mixture of spectral information from several adjacent spatial pixels of the scene. It is noted that the proposed system encodes the spectral images in both of spectral and spatial domains by using the LCTF and DMD.

In order to improve the reconstruction performance, multiple snapshots are conducted to increase the amount of measurements. For each snapshot, the coded aperture and center wavelength of LCTF are determined. In the measurement stage, we first fix the coded aperture, and then tune the LCTF to capture several snapshots in different spectral channels. Based on the sparsity assumption of the hyperspectral data, super-resolution hyperspectral images can be reconstructed from the compressive measurements using inverse optimization algorithms. Benefiting from the high resolution of the coded aperture, the spatial resolution of reconstructed hyperspectral images will be higher than that of the detector. On the other hand, the number of reconstructed spectral bands can be much more than the LCTF spectral channels since the actual transmission functions of the LCTF are used to take into account the spectral coding effects on the hyperspectral images. The two-step iterative shrinkage/thresholding (TwIST) algorithm is used to solve for the reconstruction problem, since it is computationally efficient and robust to the variation of parameters [29]. Other convex programming algorithms could also be used [30–32]. The proposed spectral imaging system and approaches are verified and assessed by a set of simulations and experiments. The results show that the proposed method can effectively improve the spectral and spatial resolutions of traditional LCTF-based spectral imaging system. Furthermore, the use of a low-resolution detector makes it possible to reduce the size and cost of the system.

Finally, we put forward a multi-channel spectral coding method to further increase the compression capacity of the system. During one integration time interval of the detector, the LCTF is randomly switched to encompass several spectral channels into one snapshot. By doing so, each measurement includes more abundant spectral information, thus the entire data cube can be reconstructed from less measurements. A set of simulations verify the feasibility of the multi-channel spectral coding method, while the experiments will be done in the future.

The remainder of this paper is organized as follows. Section 2 introduces the relevant fundamentals of CS theory. Section 3 describes the architecture and imaging model of the proposed hyperspectral imaging system. Section 4 develops the method to reconstruct the high-resolution hyperspectral images from the compressive measurements. Section 5 presents the simulation results. Section 6 describes our testbed and illustrates the experimental results. Section 7 proposes and discusses the multi-channel spectral coding method. Conclusions are provided in Section 8.

2. Fundamentals of CS theory

It is known that most of natural signals and images can be sparsely represented on some bases. Consider a one-dimensional signal f ∈ ℝ^N^×1, where N is the length of the signal, and ℝ^N^×1 represents the N × 1 real number space. Suppose f can be linearly represented on the basis Ψ ∈ ℝ^N^×^N, i. e.

f = Ψ θ,

where Ψ = [ψ₁, ψ₂, …, ψ_n], ψ_i ∈ ℝ^N^×1 denotes the i^th column of Ψ, and θ ∈ ℝ^M^×1 represents the coefficient vector. The signal f is called sparse or compressible if only S(S ≪ N) elements in θ are non-zero or have absolute values much larger than zero. Hereafter, Ψ and θ are referred to as the sparse basis and sparse coefficients of the signal f, respectively.

CS theory provides an efficient way to jointly measure and compress sparse signals [17, 18]. Consider a set of compressive measurements g ∈ ℝ^M^×1 given by

g = Φ f = Φ Ψ θ,

where Φ ∈ ℝ⁽^M ^× ^N⁾ (M ≪ n) is a random projection matrix. The conventional projection matrices include the Gaussian random matrix, Bernouli random matrix and so on. The original signal f can be recovered from g with high probability if the number of measurements satisfies

M ⩾ C \cdot μ^{2} (Φ, Ψ) \cdot S \cdot \log N,

where C is the oversampling factor,

μ (Φ, Ψ) \in [1, \sqrt{N}]

is the mutual coherence between Φ and Ψ. In particular, the mutual coherence metric is defined as

μ (Φ, Ψ) = \sqrt{N} \cdot \max_{m, n} | 〈 φ_{m}, ψ_{n} 〉 |,

where φ_m is the m^th row of Φ, and ψ_n is the n^th column of Ψ. The mutual coherence µ indicates the largest correlation between the rows of Φ and the columns of Ψ. In practice, the compressive measurements may be contaminated by noise, and then the original signal can be recovered by solving for the following optimization problem [19]

\hat{θ} = \arg \underset{θ}{m i n} {‖ θ ‖}_{1} subject to {‖ g - A θ ‖}_{2} ⩽ ε,

where ║·║₁ and ║·║₂ represent the ℓ₁-norm and ℓ₂-norm, respectively. A = ΦΨ is called the sensing matrix, and ε is the bound of noise. The optimization problem in Eq. (5) can be efficiently solved using the standard convex programming algorithms [29–31].

3. Modelling of CATF spectral imaging system

The sketch of the proposed CATF spectral imager is shown in Fig. 1(b). The spectral scene is represented by a three-dimensional data cube with two spatial dimensions along x and y axes, and one spectral dimension along λ axis. An imaging lens is used to form the image of the scene on the LCTF plane. The LCTF serves as a band-pass spectral filter over all spatial pixels of the incident scene. After the LCTF, the narrow-band spectral image is projected by a relay lens onto the coded aperture. The coded pattern is realized by the DMD that could be patterned rapidly through changing the reflection angles of the micromirrors. The narrow-band spectral image is modulated by the coded aperture in spatial domain, and then multiplexed and projected onto a low-resolution detector array via the second relay lens. We use a monochrome CMOS detector to capture the compressive measurements of the incident hyperspectral data since it has a high output speed and a small size.

In the proposed system, the resolution of the coded aperture is higher than that of the detector. Thus, the spatial resolution of the reconstructed hyperspectral images is determined by the resolution of coded aperture. On the other hand, super-resolution in the spectral domain is achieved by using an accurate transmission model of LCTF that will be described shortly. In the following, we describe the imaging model of the proposed system, where the impulse response of the imaging system is approximated by an ideal delta function.

3.1. Spectral filtering using LCTF

In the proposed system, the LCTF is used to filter the incident hyperspectral data, and collect the spectral images in different spectral channels. Conventionally, the LCTF is approximated as an ideal spectral filter with minimal bandwidth. The out-of-band transmittance of LCTF is approximated to zero, and the output of LCTF is recognized as a monochromatic spectral slice among the hyperspectral images [8]. However, the output of LCTF is indeed a multi-spectral image although the spectrum is constrained within a narrow band. The assumption of infinitesimal bandwidth will inevitably induce inaccuracy in the reconstruction stage. Thus, in this paper we adopt the actual transmission functions of the LCTF to rigorously depict its filtering characteristics.

Most current LCTFs are designed based on Lyot filters, which have the transmission functions as following [10]

T_{s} (λ) = {[\frac{\sin (2^{n} π d Δ n / λ)}{2^{n} \cdot \sin (π d Δ n / λ)}]}^{2},

where n is the number of retarders, d is the thickness of the thinnest retarder, and Δn is the birefringence. The transmission function reaches its local maxima when dΔn/λ is an integer. The birefringence of liquid crystals can be controlled by the applied voltage, which is used to adjust the center wavelength of LCTF within its free spectral range. In Fig. 2, the ideal and actual transmission functions of LCTFs are illustrated by the black dashed line and the red solid curve, respectively. The bandwidth of the ideal LCTF is considered infinitesimal, but the actual bandwidth of typical LCTF devices is in the order of 10–30 nanometers range. In this paper, the bandwidth of LCTF is defined as the full width at half maximum (FWHM) of the transmission function.

Fig. 2 The ideal and actual transmission functions of LCTFs.

Download Full Size | PDF

Figure 3 illustrates the spectral filtering characteristics of the LCTF. The LCTF is used to scan and capture the spectral images in different spectral channels by tuning the center wavelength in turn. Every spectral channel spans a specific spectral range around a particular center wavelength. Let f₀(x, y, λ) represent the spectral density of the incident scene. Suppose that we capture the spectral images from L different channels. The center wavelength of the l^th (l = 1, 2, …, L) channel is denoted by λ_l, and the corresponding transmission function is represented by $T_{s}^{l} (λ)$ . Thus, the spectral density in the l^th channel is formulated as

f_{1}^{l} (x, y, λ) = f_{0} (x, y, λ) \cdot T_{s}^{l} (λ) .

Fig. 3 The spectral filtering characteristics of the LCTF.

Download Full Size | PDF

According to Eq. (7), the LCTF can be used to modulate the hyperspectral images in the spectral domain.

3.2. Spatial coding using DMD

As shown in Fig. 1(b), the narrow-band spectral images passing through the LCTF will be modulated in the spatial domain by the coded aperture. For a given coded aperture pattern, we can obtain one snapshot of the encoded spectral images on the detector. However, a single snapshot is usually not adequate to accurately reconstruct the hyperspectral data. To solve this problem, multiple snapshots are carried out to improve the reconstruction performance as described below. Firstly, we fix the coded aperture pattern and tune the LCTF to obtain a series of snapshots in different spectral channels until all of the predefined center wavelengths are scanned. Afterwards, we load another coded aperture pattern, and repeat the aforementioned procedure. Note that in each snapshot, the coded aperture and the center wavelength of LCTF are determined, which is different from the multi-channel spectral coding method in Section 7.

Figure 4 shows the spatial coding process using K different coded aperture patterns. The transmission function of the kth coded aperture pattern is denoted by $T_{c}^{k} (x, y)$ , where k = 1, 2, …, K. In the proposed system, the coded aperture is implemented by the DMD with higher spatial resolution than the detector. The megapixel DMD generates the pre-defined binary coded aperture patterns using the individually controllable micromirrors. Each micromirror can be oriented to +12° or −12° diagonally, and reflects the incident light into different directions. Then, the encoded narrow-band spectral images by DMD is represented by

\begin{matrix} f_{2}^{l, k} (x, y, λ) = f_{1}^{l} (x, y, λ) \cdot T_{c}^{k} (x, y) \\ = f_{0} (x, y, λ) \cdot T_{s}^{l} (λ) \cdot T_{c}^{l} (x, y), \end{matrix}

where l indexes the spectral channel number (l = 1, 2, …, L), and k indexes the snapshot number (k = 1, 2, …, K).

Fig. 4 The spatial coding and compressive measurements of the spectral images.

Download Full Size | PDF

3.3. Compressive measurements on the detector

As shown in Fig. 1(b), the second relay lens is placed on one of the reflected light paths of the DMD to collect and project the encoded spectral images onto the detector plane. For the k^th snapshot in the l^th spectral channel, the light intensity on the detector plane can be obtained by integrating Eq. (8) along the wavelength:

\begin{matrix} g^{l, k} (x, y) = \int f_{2}^{l, k} (x, y, λ) d λ \\ = \int f_{0} (x, y, λ) \cdot T_{s}^{l} (λ) \cdot T_{c}^{k} (x, y) d λ . \end{matrix}

In order to facilitate the numerical analysis and computation, the imaging model is reformulated into a discrete form. All of the data are discretized in both spatial and spectral domains. Let $F_{n_{x}, n_{y}, n_{λ}}$ denote the discrete form of the incident spectral density f₀(x, y, λ), where n_x, n_y and n_λ represent the coordinates of the hyperspectral data cube in the x, y and λ axes, respectively. In particular, n_x = 1, 2, …, N_x, n_y = 1, 2, …, N_y and n_λ = 1, 2, …, N_λ. That is, the dimension of the coded aperture is N_x × N_y, and the hyperspectral data cube includes N_λ spectral bands. Let δ_c and δ_d represent the pixel sizes of the coded aperture and the detector, respectively. Assume δ_d = Rδ_c, where R is a positive integer larger than 1. Since R > 1, the spatial resolution of the proposed system is determined by the resolution of the coded aperture. The compression ratio of the spectral images in the spatial domain is γ_c = R².

Taking into account the effects of both LCTF and DMD, the overall transmission function of the proposed system can be formulated as

\begin{matrix} T^{l, k} (x, y, λ) = T_{s}^{l} (λ) \cdot T_{c}^{k} (x, y) \\ = \sum_{n_{x}, n_{y}} T_{n_{x}, n_{y}, n_{λ}}^{l, k} \cdot rect (\frac{x}{δ_{c}} - n_{x}, \frac{y}{δ_{c}} - n_{y}), \end{matrix}

where

T_{n_{x}, n_{y}, n_{λ}}^{l, k}

represents the transmittance corresponding to the voxel with coordinate (n_x, n_y, n_λ) in the hyperspectral data cube, and rect(·) is the rectangular function. The superscripts k and l represent the k^th coded aperture pattern applied in the l^th spectral channel. It is noted that the spectral resolution of the reconstructed hyperspectral images depends on the discretization precision of the LCTF transmission function

T_{s}^{l} (λ)

.

The dimension of the detector is M_x×M_y, where M_x = N_x/R and M_y = N_y/R. The compressive measurement on the (m_x, m_y)^th pixel of the detector can be calculated by integrating g^l,k(x, y) within the pixel region

\begin{matrix} G_{m_{x}, m_{y}}^{l, k} = \iint g^{l, k} (x, y) \cdot rect (\frac{x}{δ_{d}} - m_{x}, \frac{y}{δ_{d}} - m_{y}) d x f y \\ = ∭ \sum_{n_{x}, n_{y}} T_{n_{x}, n_{y}, n_{λ}}^{l, k} \cdot rect (\frac{x}{δ_{c}} - n_{x}, \frac{y}{δ_{c}} - n_{y}) \\ \cdot rect(\frac{x}{δ_{d}} - m_{x}, \frac{y}{δ_{d}} - m_{y}) \cdot f_{0} (x, y, λ) d x d y d λ \\ = \sum_{n_{x} = R (m_{x} - 1) + 1}^{m_{x} R} \sum_{n_{y} = R (m_{y} - 1) + 1}^{m_{y} R} \sum_{n_{λ} = 1}^{N_{λ}} F_{n_{x}, n_{y}, n_{λ}} \cdot T_{n_{x}, n_{y}, n_{λ}}^{l, k}, \end{matrix}

where g^l,k(x, y) is described in Eq. (9)m_x = 1, 2, …, M_x, m_y = 1, 2, …, M_y, and

F_{n_{x}, n_{y}, n_{λ}}

represents the voxel with coordinate (n_x, n_y, n_λ) in the hyperspectral data cube.

4. Reconstruction method of hyperspectral images

Assume $g \in ℝ^{K \cdot L \cdot M_{x} \cdot M_{y} \times 1}$ is the raster-scanned vector of the compressive measurements G, and $f \in ℝ^{N_{x} \cdot N_{y} \cdot N_{λ} \times 1}$ is the raster-scanned vector of the hyperspectral data cube F. The imaging model of the proposed system is reformulated as following

g = Φ f,

where

g = {[g_{1}^{T}, g_{2}^{T}, \dots, g_{M_{x} \cdot M_{y}}^{T}]}^{T}

, and g_i ∈ ℝ^K^·^L^×1 represents the measurements on the i^th detector pixel across the L spectral channels using K different coded apertures.

Φ \in ℝ^{(K \cdot L \cdot M_{x} \cdot M_{y}) \times (N_{x} \cdot N_{y} \cdot N_{λ})}

is the transmission matrix of the entire system.

f = {[f_{1}^{T}, f_{2}^{T}, \dots, f_{N_{x} \cdot N_{y}}^{T}]}^{T}

, and

f_{i} \in ℝ^{N_{λ} \times 1} (i = 1, 2, \dots, N_{x} \cdot N_{y})

denotes the spectrum of the i^th spatial pixel across the N_λ spectral bands. The transmission matrix Φ, including the effects of both LCTF and coded apertures, can be calculated as

Φ = Φ_{x y} \otimes Φ_{λ},

where ⊗ represents the Kronecker product, and

Φ_{x y} \in ℝ^{(K \cdot M_{x} \cdot M_{y}) \times (N_{x} \cdot N_{y})}

is the spatial transmission matrix of the coded apertures. The structure of the matrix Φ_xy is given by

Φ_{x y} = [\begin{matrix} Φ_{x y}^{1} & 0_{K \times R^{2}} & \dots & 0_{K \times R^{2}} \\ 0_{K \times R^{2}} & Φ_{x y}^{2} & \dots & 0_{K \times R^{2}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0_{K \times R^{2}} & 0_{K \times R^{2}} & \dots & Φ_{x y}^{M_{x} \times M_{y}} \end{matrix}],

where

0_{K \times R^{2}} \in ℝ^{K \times R^{2}}

is a zero matrix, and

Φ_{x y}^{i} \in ℝ^{K \times R^{2}}

is the spatial transmission matrix corresponding to the i^th detector pixel. In this paper, a sparse random matrix is employed as

Φ_{x y}^{i}

[33, 34]. On each column of

Φ_{x y}^{i}

, we randomly select a fix number of elements and set them to one. The sparse random matrix contains a few number of nonzero elements, so that it can reduce the computational complexity of the reconstruction algorithm. Assume

Φ_{x y}^{i} = {[{(ϕ_{x y}^{i, 1})}^{T}, {(ϕ_{x y}^{i, 2})}^{T}, \dots, {(ϕ_{x y}^{i, K})}^{T}]}^{T}

, where

ϕ_{x y}^{i, k} \in ℝ^{1 \times R^{2}}

is the k^th row of

Φ_{x y}^{i}

For the k^th snapshot,

ϕ_{x y}^{i, k}

indicates the relationship between the i^th detector pixel and the corresponding R × R subregion on the spectral images. We can stack

ϕ_{x y}^{i, k}

into a R × R matrix denoted by

Γ_{x y}^{i, k}

. It is noted that

Γ_{x y}^{i, k}

represents the R × R subregion on the k^th coded aperture. The entire coded aperture is composed of all of the submatrices

Γ_{x y}^{1, k}

,

Γ_{x y}^{2, k}, \dots, Γ_{x y}^{M_{x} \cdot M_{y}, k}

. In addition, 0 and 1 represent the block and unblock pixels, respectively. In Eq. (13),

Φ_{λ} \in ℝ^{L \times N_{λ}}

is the spectral transmission matrix of the LCTF. The l^th row of Φ_λ is generated by discretizing

T_{s}^{l} (λ)

into N_λ points. Each row represents one LCTF channel which has a corresponding Δn in Eq. (6).

Figure 5 provides an illustrative example of the transmission matrix Φ for a hyperspectral data cube with the spatial dimensions N_x = N_y = 4 and spectral dimension N_λ = 16. The number of coded apertures and spectral channels are K = 5 and L = 4, respectively. The dimensions of detector are M_x = M_y = 1. Therefore, the sizes of Φ_xy, Φ_λ and Φ are 5 × 16, 4 × 16 and 20 × 256, respectively.

Fig. 5 An example of the transmission matrix for a hyperspectral data cube with N_x = N_y = 4, N_λ = 16, K = 5, L = 4 and M_x = M_y = 1.

Download Full Size | PDF

It is noted that the dimension of the hyperspectral date cube is N_x · N_y · N_λ, and the dimension of the compressive measurements is K · L · M_x · M_y. Thus, the overall compression ratio of the proposed system is

γ = \frac{N_{x} \cdot N_{y} \cdot N_{λ}}{K \cdot L \cdot M_{x} \cdot M_{y}} = \frac{R^{2} \cdot N_{λ}}{K \cdot L} .

Since the dimension of compressive measurements is much lower than that of hyperspectral data cube, Eq. (12) is an underdetermined system of equations, which does not have the unique solution. However, based on the CS theory we can reconstruct the hyperspectral images from the compressive measurements under the sparsity assumption. It is known that the hyperspectral images can be sparsely represented on some bases, such as the discrete cosine transform (DCT) basis, various discrete wavelet transform (DWT) bases and so on. Assume that the hyperspectral data can be formulated as

f = Ψ θ,

where

Ψ \in ℝ^{(N_{x} \cdot N_{y} \cdot N_{λ}) \times (N_{x} \cdot N_{y} \cdot N_{λ})}

represents the sparse basis,

θ \in ℝ^{(N_{x} \cdot N_{y} \cdot N_{λ}) \times 1}

is the corresponding sparse coefficient vector. The sparse basis Ψ is defined as Ψ = Ψ₁ ⊗ Ψ₂ ⊗ Ψ₃, where

Ψ_{1} \otimes Ψ_{2} \in ℝ^{(N_{x} \cdot N_{y}) \times (N_{x} \cdot N_{y})}

is the two-dimensional Haar DWT basis, and

Ψ_{3} \in ℝ^{N_{λ} \times N_{λ}}

is the one-dimensional DCT basis [35]. Substituting Eq. (16) into Eq. (12) and taking into account the measurement noise on the detector, we have

g = Φ Ψ θ + ω = A θ + ω,

where A = ΦΨ is referred to as the sensing matrix, and

ω \in ℝ^{K \cdot L \cdot M_{x} \cdot M_{y} \times 1}

is the noise vector. Thus, the underlying hyperspectral images can be reconstructed by solving for the following l₁-norm minimization problem:

\hat{θ} = \arg \min_{θ} {‖ θ ‖}_{1} subject to {‖ g - A θ ‖}_{2} ⩽ ε,

where ε is the bound of the noise. In this paper, we use the TwIST algorithm to solve the reconstruction problem in Eq. (18), since it is computational efficient and robust to the variation of parameters [29]. Other optimization algorithms could be used as well [30–32]. The algorithms are implemented by MATLAB R2016a. All of the computations are carried out on a server with Intel Xeon E3-1505M v5 2.8GHz processor, and 64GB memory.

5. Simulation results

This section verifies the proposed method by a set of simulations using the real hyperspectral data. The original hyperspectral images are extracted from the hyperspectral database of Ben-Gurion University Interdisciplinary Computational Vision Laboratory [36]. Then, we truncate the hyperspectral images from the database to reduce the data volume and to fit the spectral sensitivity range of the detector. The truncated hyperspectral images consist of 400 × 400 pixels in the spatial domain, and 170 spectral bands from 500nm to 710nm. The spectral width of each band is about 1.24nm. In addition, we normalize the intensity of the hyperspectral images into the range of [0, 1]. The transmission function in Eq. (6) is used to represent the spectral filtering characteristics of the LCTF (as shown in Fig. 6).

Fig. 6 The LCTF transmission functions used in the simulations. The center wavelengths of LCTF are ranged from 500nm to 710nm with the interval of 10nm.

Download Full Size | PDF

The top row of Fig. 7 shows the original spectral images at 510nm, 600nm and 690nm. The middle row of Fig. 7 shows the reconstructed spectral images using the proposed CS method based on the CATF spectral imager in Fig. 1(b). The LCTF is used to filter the spectral images into 22 spectral channels. The center wavelengths of the spectral channels are ranged from 500nm to 710nm with the interval of 10nm. The spectral compression ratio of the system is γ_s = N_λ/L = 170/22 ≈ 7.73. The dimension of the coded aperture is 400 × 400, and the dimension of the detector is 50 × 50. Thus, the spatial compression ratio of the system is γ_c = R² = (400/50)² = 64. Since we capture K = 25 snapshots in each spectral channel, the overall compression ratio of the system is γ = γ_s · γ_c/K = 7.73·64/25 ≈ 19.79. The measurement noise on the detector is white Gaussian noise with signal to noise ratio (SNR) of 40dB. The reconstructed spectral images consist of 170 spectral bands and 400 × 400 spatial pixels. The peak signal to noise ratios (PSNR) of the reconstructed spectral images are presented in the middle row of Fig. 7, and the average PSNR over all spectral bands is 28.01dB. It is shown that the proposed methods can recover the high-resolution spectral images with promising image quality.

Fig. 7 (Top) the original spectral images, (middle) the reconstructed spectral images using the CATF system, and (bottom) the reconstructed spectral images using the traditional system.

Download Full Size | PDF

As a comparison, the bottom row of Fig. 7 illustrates the spectral images acquired by the traditional LCTF-based system in Fig. 1(a). Similar to the simulations in the middle row of Fig. 7, the LCTF is used to extract the spectral images from 22 spectral channels, and the dimension of the detector is 50 × 50. The traditional system falls short to improve the spectral resolution since it approximates the LCTF as an ideal spectral filter with minimal bandwidth. In addition, the traditional system obtains the spectral images using the low-resolution detector instead of enhancing the resolution using CS methods. In order to make a comparison with the proposed system, we extend the 22 spectral bands to 170 spectral bands using the linear interpolation method. Thus, the spectral images acquired by the traditional system consist of 50 × 50 spatial pixels and 170 spectral bands. The PSNRs of the reconstructed spectral images are presented in the bottom row of Fig. 7, and the average PSNR is 21.69dB. Compared to the traditional method, the proposed method can improve PSNR by 6.32dB on average.

Figure 8 compares the reconstructed spectra of three representative pixels on the scene using different systems. Figure 8(a) shows the RGB image of the scene. The three representative pixels are located at P1, P2 and P3, respectively. These three pixels come from different colorful regions, thus having different spectral characteristics. Figures 8(b)–8(d) present the spectra on P1, P2 and P3, respectively. In particular, the black solid lines represent the original spectra. The blue dash-dotted lines and red dashed lines represent the reconstructed spectra using the traditional and proposed methods, respectively. For the proposed method, the PSNRs of the reconstructed spectra on P1, P2 and P3 are 28.47dB, 30.27dB and 31.35dB, respectively. On the other hand, the PSNRs of the reconstructed spectra obtained by the traditional method are 17.17dB, 19.15dB and 26.55dB, respectively. It is observed that the proposed method leads to superior reconstruction performance over the traditional method.

Fig. 8 The original and reconstructed spectra of three representative pixels on the scene using different methods. (a) The RGB image of the scene, and the spectra on the pixels (b) P1, (c) P2 and (d) P3.

Download Full Size | PDF

6. Experimental results

This section describes the testbed of CATF spectral imager developed by our group. Figure 9 shows the picture of the testbed, which consists of an imaging lens, an LCTF, two relay lenses, a DMD and a monochromatic CMOS camera. The object is illuminated by a broadband annular light source. The light rays emitted from the object will successively transmit through the imaging lens, the LCTF, the first relay lens and then projected onto the DMD. After that, the second relay lens is used to collect the spectral images onto the CMOS camera. The focal lengths of imaging lens and relay lenses are 50mm and 900mm, respectively. The wavelength range of the LCTF is 500nm–710nm. At the center wavelength of 500nm, the bandwidth of LCTF is 10nm. At the center wavelength of 710nm, the bandwidth of LCTF is 15nm. The DMD (DLP9500) consists of 1920 × 1080 pixels with pixel pitch of 10.8µm. The resolution of CMOS is 1024 × 1024, and the pixel pitch is 11.0µm.

Fig. 9 Testbed of the proposed CATF spectral imager.

Download Full Size | PDF

In the experiment, the calibration process is necessary to compensate for the non-ideal transmittance of the coded aperture, the spectral response of the optical elements, and the detector response. For all of the coded aperture patterns used in the experiment, we measured their patterns on the detector when the object is replaced by a diffuse plate. Then, we acquired the spectral response of each component in the system by a grating spectrometer (HR4000 Ocean Optics), including the spectral transmittance of each lens, the transmission function of the LCTF under each spectral channel, and the spectral reflectance of the DMD as all of its pixels turned on. In addition, the spectral response of the detector was obtained from its specification. In the reconstruction process, the normalized detected images of the coded apertures were used as the actual spatial transmission matrix. Then, the spectral transmittance or reflectance of all components and the detector response were included in the spectral transmission matrix of the system.

In the measurement stage, 25 different coded aperture patterns are used to conduct the spatial modulation. For each coded aperture pattern, the LCTF is switched for 22 times to capture the snapshots in 22 spectral channels. The center wavelengths of the spectral channels are from 500nm to 710nm with an interval of 10nm. On the DMD, we choose an area including 400 × 400 pixels to implement the coded aperture. On the detector, we merge 8 × 8 pixels into one macro-pixel. The ratio of pixel pitches between the detector and the coded aperture is 8, that is δ_d = 8δ_c. Thus, we use 50 × 50 macro-pixels on the detector to collect the measurements. The reconstructed spectral images consist of 170 spectral bands and 400 × 400 spatial pixels.

As a comparison, we also construct the testbed of the traditional spectral imaging system illustrated by Fig. 1(a). The structure of the traditional system is similar to the proposed system. However, in the traditional system all pixels on the DMD will be turned on. Thus, the DMD will directly reflect the incident light onto the second relay lens without spatial coding. The spectral images obtained by the traditional system consist of 22 spectral bands and 50 × 50 spatial pixels. After that, we extend the 22 spectral bands to 170 spectral bands using linear interpolation method.

Figure 10 shows the actual spectral transmission functions of the entire imaging systems corresponding to different center wavelengths. The actual transmission functions are measured by a grating spectrometer, and the discretization precision of the transmission functions is set to 1.24nm. The actual imaging system is not ideal, and its spectral transmittance is affected by the LCTF and other components, such as the lenses, the DMD, and so on.

Fig. 10 The actual spectral transmission functions of the entire system corresponding to different center wavelengths.

Download Full Size | PDF

Figure 11 shows the reconstructed spectral images at 530nm, 580nm and 630nm. The reconstructed images using the proposed system contain 400 × 400 pixels in the spatial domain, while the reconstructed images using the traditional system only contain 50 × 50 pixels. That is because the traditional system does not use the CS method to enhance its spatial resolution.

Fig. 11 The reconstructed spectral images using (top) the CATF system and (bottom) the traditional system.

Download Full Size | PDF

Figure 12 shows the original and reconstructed spectra corresponding to three representative pixels on the scene, where the representative pixels are located at P1, P2 and P3 in Fig. 12(a). The original spectra are measured by the grating spectrometer. Figures 12(b)–12(d) present the reconstructed spectra. For the proposed method, the PSNRs of the reconstructed spectra on P1, P2 and P3 are 31.43dB, 26.48dB and 30.48dB, respectively. For the traditional method, the PSNRs of the reconstructed spectra are 22.85dB, 20.14dB and 22.74dB, respectively. The experimental results show that the proposed method can effectively improve the reconstruction performance in contrast to the traditional method. The CATF system could be exploited in the situations, where the low-resolution detector is available to acquire and reconstruct the high-resolution spectral images.

Fig. 12 The spectra of three representative pixels on the scene obtained by the testbed. (a) The RGB image of the scene, and the spectra on the pixels of (b) P1, (c) P2 and (d) P3.

Download Full Size | PDF

It is worth noting that there are several factors influencing the imaging performance of our current experimental testbed, such that the simulation results look better than the actual experimental results. First, the LCTF and DMD in the system will attenuate the incident light energy. The low utilization rate of light energy is an inherent defect of the narrow-band filter device, and the actual peak transmittance of the LCTF under each spectral channel is 10%–20%. Besides, the transmission rate of DMD is designed to about 50% for the random coded apertures. To achieve higher image intensity levels, the integration time of the detector has to be increased, which can result in a lower signal to noise ratio attributed to the noise of the detector. Second, the system is composed by several separated elements without encapsulation. The DMD and the detector will be affected by the stray light, thus degrading the reconstruction performance. Finally, the impulse response of the system is considered an ideal delta function, which is not the actual case. The impulse response of the actual experiment system will blur the detected images and increase reconstruction errors.

7. Multi-channel spectral coding method

In the CATF spectral imager mentioned above, the center wavelength of LCTF is fixed during one snapshot. Thus, only a narrow-band spectral channel is involved in one snapshot. The high reconstruction quality requires a number of snapshots to provide enough spectral information at the cost of reducing the compression capacity of the system, where the compression ratio is defined in Eq. (15). This section proposes a multi-channel spectral coding method to further improve the compression capacity and maintain the desirable reconstruction performance at the same time. During one integration time interval of the detector, we fix the coded aperture pattern, but randomly switch the LCTF for several times to collect the information from multiple spectral channels in one snapshot. Therefore, each measurement includes more abundant spectral information, thus the spectral images can be reconstructed from less measurements. In order to distinguish, we refer to the spectral coding method in Section 3.1 as the single-channel spectral coding method.

Suppose that Q snapshots are captured under every given coded aperture pattern. In the qth (q = 1, 2, …, Q) snapshot, the center wavelength of LCTF is randomly switched by P_q times to collect the information from P_q different spectral channels. Thus, the spectral transmission function of the qth (q = 1, 2, …, Q) snapshot can be written as

T_{P_{q}} (λ) = \sum_{i = 1}^{P_{q}} {[\frac{\sin (2^{n} π d Δ n_{i} / λ)}{2^{n} \cdot \sin (π d Δ n_{i} / λ)}]}^{2},

where different Δn_i corresponds to different spectral channel. Compared to the single-channel spectral coding in Eq. (6), Eq. (19) encompass the transmission functions from different spectral channels. To make a fair comparison between the two spectral coding methods, the total number of involved spectral channels should be the same. That is,

\sum_{q = 1}^{Q} P_{q} = L

, where L is described above Eq. (7). The spectral transmission matrix of the system is denoted as

Φ_{λ} \in ℝ^{Q \times N_{λ}}

, where the q^th row of Φ_λ is generated by discretizing

T_{P_{q}} (λ)

into N_λ points. Figure 13 shows an example of the transmission matrix Φ under the multi-channel spectral coding framework. The dimensions of the hyperspectral data cube are N_x = N_y = 4 and N_λ = 16. The dimensions of the detector are M_x = M_y = 1. The number of coded aperture patterns used is K = 5. Under each coded aperture pattern, we make Q = 2 snapshots by switching the LCTF. Compared to Fig. 5, the number of rows in Φ reduces by half since the compression ratio in spectral dimension has doubled.

Fig. 13 An example of the transmission matrix using multi-channel spectral coding method.

Download Full Size | PDF

Next, we verify the proposed multi-channel spectral coding method by simultions, where the same hyperspectral images are used as in Section 5. All of the parameters are the same as those in Fig. 7(b) except for the spectral transmission matrix Φ_λ. The center wavelengths of the LCTF are ranged from 500nm to 710nm but multiple spectral channels are collected during each snapshots. In particular, we randomly divide the 22 spectral channels into Q groups, and every group includes the same number or almost the same number of spectral channels. In the measurement stage, each snapshot collects information from one group of spectral channels, thus each spectral channel is collected only once. For example, if Q = 11, then P_q = 2 for all q. If Q = 8, then P_q = 2 for q ≤ 2 and P_q = 3 for q > 2. Table 1 shows the maximum, minimum and average PSNRs of the reconstructed spectral images for different number of Q. For each number of Q, we calculate the PSNRs on average by repeating the simulations for 15 times.

Table 1. The maximum, minimum and average PSNRs of the reconstructed spectral images for different number of Q.

View Table

It is observed that the PSNRs decrease along with the decrement of Q since the compression ratio of the system increases. On the other hand, the compression ratio is decreased as the Q increases, which benefits the improvement of PSNRs. Thus, Table 1 shows the intrinsic tradeoff between the compression ratio and the reconstruction performance. It is noted that the proposed multi-channel spectral coding method may improve the compression capacity of the system, while maintain the quality of the reconstructed spectral images. Inspired by these simulations, it is possible to obtain better reconstruction results by optimizing the multi-channel spectral coding framework, where the optimal spectral transmission matrix can be solved based on the restricted isometry property in CS theory [26, 37]. Due to the length limitation of the paper, the optimization of spectral coding and the experiments will be done in the future.

8. Conclusion

This paper proposed a super-resolution hyperspectral imaging system based on the compressive sensing framework. The LCTF is used to filter and modulate the spectra of the scene. In addition, the coded aperture is used to modulate the spectral images in the spatial domain. Based on the sparsity assumption, the super-resolution hyperspectral images can be reconstructed from a small set of compressive measurements collected by the low-resolution detector. The testbed was established to verify the proposed spectral imaging method. The simulations and experiments show that the proposed method can achieve superior reconstruction performance of spectral images compared to the traditional method. Finally, the multi-channel spectral coding method was proposed and the future work was discussed.

Funding

National Natural Science Foundation of China (NSFC) (61527802, 61371132, 61471043); National Science Foundation (NSF) (VEC 1538950).

References

1. H. Kurosaki, H. Koshiishi, T. Suzuki, and K. Tsuchiya, “Development of tunable imaging spectro-polarimeter for remote sensing,” Adv. Space Res. 32(11), 2141–2146 (2003). [CrossRef]

2. S. C. Gebhart, R. C. Thompson, and A. Mahadevan-Jansen, “Liquid-crystal tunable filter spectral imaging for brain tumor demarcation,” Appl. Opt. 46(10), 1896–1910 (2007). [CrossRef] [PubMed]

3. S. J. Woltman, G. D. Jay, and G. P. Crawford, “Liquid-crystal materials find a new order in biomedical applications,” Nat. Mater. 6(12), 929–938 (2007). [CrossRef] [PubMed]

4. J. Beeckman, K. Neyts, and P. J. M. Vanbrabant, “Liquid-crystal photonic applications,” Opt. Eng. 50(8), 081202 (2011). [CrossRef]

5. W. Wang, C. Li, E. W. Tollner, G. C. Rains, and R. D. Gitaitis, “A liquid crystal tunable filter based shortwave infrared spectral imaging system: Calibration and characterization,” Comput. Electron. Agr. 80, 135–144 (2012). [CrossRef]

6. A. Mian and R. Hartley, “Hyperspectral video restoration using optical flow and sparse coding,” Opt. Express 20(10), 10658–10673 (2012). [CrossRef] [PubMed]

7. D. Jiang, Y. Liu, W. Cao, and Z. Guo, “A compact novel microwave liquid crystal tunable band-pass filter,” Optik 127(20), 8374–8378 (2016). [CrossRef]

8. M. A. López-Álvarez, J. Hernández-Andrés, and J. Romero, “Developing an optimum computer-designed multispectral system comprising a monochrome ccd camera and a liquid-crystal tunable filter,” Appl. Opt. 47(24), 4381–4390 (2008). [CrossRef] [PubMed]

9. G. Shabtay, E. Eidinger, Z. Zalevsky, D. Mendlovic, and E. Marom, “Tunable birefringent filters – optimal iterative design,” Opt. Express 10(26), 1534–1541 (2002). [CrossRef] [PubMed]

10. O. Aharon and I. Abdulhalim, “Liquid crystal lyot tunable filter with extended free spectral range,” Opt. Express 17(14), 11426–11433 (2009). [CrossRef] [PubMed]

11. M. Abuleil and I. Abdulhalim, “Narrowband multispectral liquid crystal tunable filter,” Opt. Lett. 41(9), 1957–1960 (2016). [CrossRef] [PubMed]

12. A. Jullien, R. Pascal, U. Bortolozzo, N. Forget, and S. Residori, “High-resolution hyperspectral imaging with cascaded liquid crystal cells,” Optica 4(4), 400–405 (2017). [CrossRef]

13. R. M. Willett, R. F. Marcia, and J. M. Nichols, “Compressed sensing for practical optical imaging systems: a tutorial,” Opt. Eng. 50(7), 072601 (2011). [CrossRef]

14. G. P. Abousleman, M. W. Marcellin, and B. R. Hunt, “Compression of hyperspectral imagery using the 3-D DCT and hybrid DPCM/DCT,” IEEE Trans. Geosci. Remote Sens. 33(1), 26–34 (1995). [CrossRef]

15. F. Rizzo, B. Carpentieri, G. Motta, and J. A. Storer, “Low-complexity lossless compression of hyperspectral imagery via linear prediction,” IEEE Signal Process. Lett. 12(2), 138–141 (2005). [CrossRef]

16. Q. Du and J. E. Fowler, “Hyperspectral image compression using jpeg2000 and principal component analysis,” IEEE Geosci. Remote Sens. Lett. 4(2), 201–205 (2007). [CrossRef]

17. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]

18. E. Candès and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse Probl. 23, 969–985 (2007). [CrossRef]

19. E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25(2), 21–30 (2008). [CrossRef]

20. 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(21), 14013–14027 (2007). [CrossRef] [PubMed]

21. A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,” Appl. Opt. 47(10), B44–B51 (2008). [CrossRef] [PubMed]

22. H. Rueda, H. Arguello, and G. R. Arce, “DMD-based implementation of patterned optical filter arrays for compressive spectral imaging,” J. Opt. Soc. Am. A 32(1), 80–89 (2015). [CrossRef]

23. H. Rueda, H. Arguello, and G. R. Arce, “Compressive spectral testbed imaging system based on thin-film color-patterned filter arrays,” Appl. Opt. 55(33), 9584–9593 (2016). [CrossRef] [PubMed]

24. T. Sun and K. Kelly, “Compressive sensing hyperspectral imager,” in Frontiers in Optics 2009/Laser Science XXV/Fall 2009 OSA Optics & Photonics Technical Digest, OSA Technical Digest (CD) (Optical Society of America, 2009), paper CTuA5.

25. Y. August, C. Vachman, Y. Rivenson, and A. Stern, “Compressive hyperspectral imaging by random separable projections in both the spatial and the spectral domains,” Appl. Opt. 52(10), D46–D54 (2013). [CrossRef] [PubMed]

26. G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S. Kittle, “Compressive coded aperture spectral imaging: an introduction,” IEEE Signal Process. Mag. 31(1), 105–115 (2014). [CrossRef]

27. I. August, Y. Oiknine, M. AbuLeil, I. Abdulhalim, and A. Stern, “Miniature compressive ultraspectral imaging system utilizing a single liquid crystal phase retarder,” Sci. Rep. 6, 23524 (2016). [CrossRef]

28. Y. Oiknine, I. August, and A. Stern, “Along-track scanning using a liquid crystal compressive hyperspectral imager,” Opt. Express 24(8), 8446–8457 (2016). [CrossRef] [PubMed]

29. J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Trans. Image Process. 16(12), 2992–3004 (2007). [CrossRef] [PubMed]

30. M. A. T. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems,” IEEE J. Sel. Topics Signal Process. 1(4), 586–597 (2007). [CrossRef]

31. S. J. Wright, R. D. Nowak, and M. A. T. Figueiredo, “Sparse reconstruction by separable approximation,” IEEE Trans. Signal Process. 57(7), 2479–2493 (2009). [CrossRef]

32. J. Tan, Y. Ma, H. Rueda, D. Baron, and G. R Arce, “Compressive hyperspectral imaging via approximate message passing,” IEEE J. Sel. Topics Signal Process. 10(2), 389–401 (2016). [CrossRef]

33. A. Gilbert and P. Indyk, “Sparse recovery using sparse matrices,” Proc. IEEE 98(6), 937–947 (2010). [CrossRef]

34. B. Bah and J. Tanner, “Vanishingly sparse matrices and expander graphs, with application to compressed sensing,” IEEE Trans. Inf. Theory 59(11), 7491–7508 (2013). [CrossRef]

35. M. F. Duarte and R. G. Baraniuk, “Kronecker compressive sensing,” IEEE Trans. Image Process. 21(2), 494–504 (2012). [CrossRef]

36. B. Arad and O. Ben-Shahar, “Sparse Recovery of Hyperspectral Signal from Natural RGB Images,” in Proceedings of Computer Vision - ECCV 2016: 14th European Conference, (Springer, 2016), pp. 19–34.

37. Y. Eldar and G. Kutyniok, Compressed Sensing: Theory and Applications (Cambridge University, 2012). [CrossRef]

Q	max/dB	min/dB	average/dB
5	24.17	21.13	22.55
6	24.55	21.28	22.75
8	24.94	22.58	23.61
11	25.81	23.54	25.04

Compressive spectral imaging system based on liquid crystal tunable filter

Abstract

1. Introduction

2. Fundamentals of CS theory

3. Modelling of CATF spectral imaging system

3.1. Spectral filtering using LCTF

3.2. Spatial coding using DMD

3.3. Compressive measurements on the detector

4. Reconstruction method of hyperspectral images

5. Simulation results

6. Experimental results

7. Multi-channel spectral coding method

8. Conclusion

Funding

References

Cited By

Figures (13)

Tables (1)

Equations (19)

Optics Express