Open Access
Add to BibSonomy
Abstract
In recent years, the demand for hyperspectral imaging devices has grown significantly, driven by their ability of capturing high-resolution spectral information. Among the several possible optical designs for acquiring hyperspectral images, there is a growing interest in interferometric spectral imaging systems based on division of aperture. These systems have the advantage of capturing snapshot acquisitions while maintaining a compact design. However, they require a careful calibration to operate properly. In this work, we present the interferometer response characterization algorithm (IRCA), a robust three-step procedure designed to characterize the transmittance response of multi-aperture imaging spectrometers based on the interferometry of Fabry-Perot. Additionally, we propose a formulation of the image formation model for such devices suitable to estimate the parameters of interest by considering the model under various regimes of finesse. The proposed algorithm processes the image output obtained from a set of monochromatic light sources and refines the results using nonlinear regression after an ad-hoc initialization. Through experimental analysis conducted on four different prototypes from the Image SPectrometer On Chip (ImSPOC) family, we validate the performance of our approach for characterization. The associated source code for this paper is available from Zenodo (http://dx.doi.org/10.5281/zenodo.7978514).
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The demand for imaging spectrometers, also known as hyperspectral (HS) cameras, has experienced significant growth in recent years. This surge in popularity can be attributed to their outstanding ability to capture high-resolution spectral information, especially in comparison to classic multispectral devices. These cameras find applications in various fields, such as astronomy, precision agriculture, molecular biology, biomedical imaging, geosciences, physics, and surveillance [1–4]. Of particular importance is their role in accurately measuring atmospheric gases, which is vital for climate change monitoring, air quality studies, and compliance to regulatory requirements [5].
Traditional imaging spectrometers that rely on scanning mechanism, such as whiskbroom and pushbroom, face limitations in capturing spatially varying scenes and are forced to make compromises between spectral and spatial resolution [6]. Consequently, significant research efforts have been recently dedicated to the development and production of computational spectral imaging systems. These systems aim to enhance spectral, spatial, and temporal resolution and operate by encoding hyperspectral information in low-dimensional projected domains. However the retrieval of the full spectral and spatial HS datacube requires the application of sophisticated reconstruction algorithms [7,8].
In this paper, we focus our attention on the characterization of the transmittance response of multi-aperture interferometric imaging spectrometers [10]. This class of instruments includes miniaturized snapshot acquisition systems for HS imagery, whose optical design consists of a matrix of microlenses and a staircase-shaped optical plate superposed to a focal plane array. Figure 1 shows the optical design of one of such devices, known as Image SPectrometer On Chip (ImSPOC) [5,11–14].
Fig. 1. Optical concept of multi-aperture interferometric imaging spectrometers and examples of ImSPOC devices. (a) Cross-section view of the optical concept for an ImSPOC imaging system. In the pictured design, the interferometers are air cavities of different thickness carved within a glass optical plate and coated with a reflective layer in titanium dioxide (TiO2). An external spectral filter can be also used to limit the input wavenumber range of the device for certain applications, such as NO2 gas detection [9]. (b), (c) Examples of Image SPectrometer On Chip (ImSPOC) prototypes; their specific characteristics are described in the experimental section.
Figure 2 illustrates an example of acquisition. The acquired image is composed of several subimages, with each subimage being the result of filtering the incident radiance using the transmittance response of a unique interferometer/microlens unit. The set of readings obtained from the various subimages, arranged in ascending order of interferometer thicknesses, can be viewed as a sampled representation of a continuous interferogram linked to the spectrum of the particular region of the scene being viewed by the device.
In comparison to traditional hyperspectral cameras, multi-aperture interferometric imaging spectrometers offer several advantages such as snapshot acquisitions, compact dimensions, while preserving competitive spectral resolutions ranging from 5 to 10 nm. However, they do face limitations in terms of spatial resolution and potential parallax effects.
Fig. 2. Visual representation of the operating principle of a multi-aperture interferometric imaging spectrometer. In this example, the acquisition is a $512 \times 640$ image taken by an ImSPOC device. The image is filtered by 70 interferometer/microlens units, resulting in the scene being split into multiple subimages. Matching pixels across these subimages are then arranged to form an interferogram for each ground point, which is subsequently processed to reconstruct the spectrum [15,16]. The purpose of this figure is purely illustrative.
The identification of the image formation model for these devices is a crucial step that boils down to the estimation of the instrument optical transformation. Furthermore, regular calibration of the instrument becomes essential to maintain up-to-date device characterizations. This is especially relevant when considering potential changes in the instrument’s physical properties over time, resulting from factors like instrument aging or variations in acquisition conditions, such as temperature. As an example scenario, due to the limited accuracy in either manufacturing or assembling the various device parts, the real thickness of the interferometers may be different with respect to the value they were designed for. As shown in Fig. 3, if this information is not taken into account, the interferogram samples are then placed incorrectly in domain of the optical path difference (OPD), which may cause inaccuracies on the quality of the reconstructed spectrum.
Fig. 3. Illustration of the effect of the errors in the estimation of the interferometer thickness. If the thicknesses of the optical components are known with an error (depicted by an orange dashed outline), the interferometer may sample the interferogram at a position that does not accurately correspond to the actual OPD value. In the specific case of the presented interferogram, the samples shift from the blue to the orange positions.
To address this issue, we propose in this work a general procedure for the parametric characterization of the instrument. The procedure is divided into a measurement session, where the device is illuminated with a set of flat field monochromatic sources, and an algorithmic estimation of the parameters of interest for the transmittance response, which we coin as interferometer response characterization algorithm (IRCA).
The image formation model for the devices that we aim to characterize is also recalled in this work. However, with respect to the typical formulation of the literature, we derive the formation model in terms of the parameters of interest. We formalize the mathematical model of the transmittance response of the device, specifically in terms of the OPD, reflectivity, gain, and phase shift. We also express this transmittance response under different finesse regimes. Each regime corresponds to a specific number of emerging waves in the Fabry-Perot (FP) cavity, arranged in descending order of optical power. Previous studies [16,17] implicitly characterized such devices under the assumption of 2 emerging waves, prioritizing conceptual simplicity over a precise parameter estimation. Our formulation enables us to describe previous techniques within our proposed framework, allowing for the application of similar trade-offs if desired.
Other than for ImSPOC devices, the IRCA can also be potentially employed to characterize and regularly update the calibration of various devices that exhibit a response based on the interferometry of Fabry-Perot. This includes compressive imagers [10,18] and hyperspectral imaging systems with dielectric mirrors [19,20], among others.
The IRCA is defined by a three step procedure: the overall optical gain is firstly addressed discarding any interferometric effect, then a first rough assessment of the remaining parameters is performed by casting the problem as a maximum likelihood (ML) estimation of the characteristics of a sinusoidal signal. We refine this estimation by casting the problem as a nonlinear regression and solving it with the Levenberg-Marquardt (LM) algorithm [21]. The nonlinear regression approach was also employed in other works [22], but we focus our attention here on a robust solution for optical devices whose sensors are particularly sensitive to noise, as the different parameters are made separable by imposing that their polynomial expression in terms of wavelength has a limited degree.
To summarize, the novel contributions of this work are:
- 1. The formalization of the image formation principle of multi-aperture Fabry-Perot imaging spectrometers (interferometers, lenslet, etc.).
We define within a single framework the dependency on its characteristic parameters (OPD, gain, reflectivity, phase shift) and the regimes of finesse associated to different amounts of transmitted waves;
- 2. The development of the IRCA, a procedure for the estimation of parameters for transmittance responses of devices operating as FP interferometers;
- 3. The definition of an experimental procedure for the characterization of multi-aperture interferometric imaging spectrometers, using monochromatic sources. We test the effectiveness of the proposed method on real acquisitions from four ImSPOC prototypes with different characteristics.
The article is organized as follows: in Section 2 we describe the image formation model of the multi-aperture interferometric imaging spectrometers, in Section 3 we describe the proposed spectral characterization setup and estimation algorithm, and in Section 4 we evaluate its performances and discuss its results in relation to the physics of the devices.
2. Image formation for multi-aperture interferometric imaging spectrometers
In this section, we describe the image formation model of a multi-aperture interferometric imaging spectrometer. We begin by deriving the expression of the transmittance response for a single FP interferometer/microlens unit in Section 2.1. We then specify it within our framework for different regimes of finesse in Section 2.2. Finally, we identify the parameters of interest for their characterization with their respective model in Section 2.3.
Following the literature of Fourier transform spectrometers, in this work the spectra and transmittance responses are expressed in terms of wavenumbers $\sigma$, that is as the reciprocal of the wavelengths (e.g., a wavelength of 500 nm corresponds to $\sigma =2\times 10^4$ cm−1), but the relevant plots include both wavenumber and wavelength scales. Furthermore, the vertical ordinates are appropriately labeled as normalized intensity when the intensity is scaled by its maximum value, and as mean scaled intensity when scaled by its mean value. In situations involving multiple plots, all plots are consistently scaled using the mean of the reference. For the reader’s convenience, the variables used in this paper are shown in Table 1, separated into variables for the continuous image formation model, for its parameters, for the acquisition vectors and the vector sizes. These variables will be formerly introduced when relevant to the discussion.
Table 1. Selection of variables used in this paper, grouped in their respective categories.
2.1 Optical transfer model
We want to define here the expression of the sensors’ readout in terms of the incident radiance. To this purpose, we analyze the light ray propagation within a single interferometer/lens unit of the optical system, as shown in Fig. 4(a).
Fig. 4. Light ray propagation within a multi-aperture FP imaging spectrometer. (a) Zoomed-in visualization of a single interferometer/microlens unit, showcasing the focusing effect of the light rays to the focal plane. (b) Detailed view of the Fabry-Perot cavity. The light propagation gives rise to an OPD, which corresponds to the difference between the optical paths highlighted in red and green.
By considering a scene at the optical infinity, there is a bijective correspondence between the direction of incidence $\boldsymbol{\omega }$ of the incident light and the position $\mathbf {s}$ on the focal plane. Consequently, we can express the spectral radiance $\mathcal {L}$ of the incident light either as $\mathcal {L}(\sigma,\,\boldsymbol{\omega })$ or $\mathcal {L}(\sigma,\,\mathbf {s})$. In this scenario, the $k$-th interferometer acts as a spectral filter and introduces an attenuation $\mathcal {T}_k$ of the radiant flux which varies only with the angle of incidence $\boldsymbol{\omega }$ and the wavenumber $\sigma$. As in the previous case, this can also be expressed interchangeably as $\mathcal {T}_k(\sigma,\,\boldsymbol{\omega })$ or $\mathcal {T}_k(\sigma,\,\mathbf {s})$.
Assuming no crosstalk in the formation of each subimage, the spectral flux $\Phi _{jk}(\sigma )$ received by the $j$-th sensor (i.e., a photodetector) is only due to incident light within a given $k$-th interferometer. Its expression at the focal plane is given by:
Considering that the etendue is conserved across the object and the image space, we can rewrite Eq. (1) at the input plane as:
Finally, we model the intensity level $x_{jk}$ captured by the photodetector as:
2.2 Fabry-Perot regimes of finesse
We now focus our attention on expanding the term $\mathcal {T}_k(\sigma,\,\boldsymbol{\omega })$ from Eq. (2). Let us consider a monochromatic plane wave with complex amplitude $\mathbf {E}^{[\mathrm {i}]}(\sigma )$ incident to the FP interferometer, forming an angle $\theta ^{[\mathrm {i}]}$ with the normal to the incident plane. The complex amplitude $\mathbf {E}^{[\mathrm {o}]}$ of the transmitted light can be seen as a sum of $W\rightarrow \infty$ successive emerging waves $\{\mathbf {E}_m\}_{ {m} \in [{0}, \ldots , {W-1}]}$.
Each emerging wave introduces a fixed round trip phase difference:
where $\varphi _0^{}$ defines a constant phase shift and $\delta$ defines the OPD between two consecutive emerging waves. By referring to the geometry shown in Fig. 4(b), the OPD is determined as the difference between the optical paths for a round trip inner reflection and a direct transmission. By making use of Snell’s law, simple geometrical manipulations yield:In the following, we denote as $\mathcal {T}_{k}^{[W]}(\sigma, \boldsymbol{\omega })$ the expression of $\mathcal {T}_k(\sigma, \boldsymbol{\omega })$ specific to the a generic integer amount $W$ of emerging waves. Its expression is defined by the ratio between the output and input irradiance and evaluates as follows:
For our purposes, it is also convenient to derive the mean scaled expression $\overline {\mathcal {T}}_{k}^{[W]}$ of $\mathcal {T}_{k}^{[W]}$ as:
Figure 5(a) presents the plot of the expression of $\overline {\mathcal {T}}_{k}^{[\infty ]}$ for different values of reflectivity. This visualization assumes that $\mathcal {R}$ and $n$ remain constant regardless of the wavenumbers $\sigma$. However, this assumption is rarely verified in more realistic scenarios, where variations with respect to $\sigma$ are often observed. The spacing between the peaks of the Airy distribution decreases as the thickness $d_k$ of the interferometer increases. This principle is exploited by multi-aperture devices to create different transmittance responses for different subimages as previously shown in Fig. 2.
Fig. 5. Regimes of reflectivity. (a) The figure illustrates the mean scaled theoretical transmittance response for various values of the surface reflectivity $\mathcal {R}$. (b) Depending on the desired level of accuracy, the user can establish a root mean square error threshold to identify the maximum reflectivity value (dashed line) in which the $W$-wave model of Eq. (7b) exhibits behavior indistinguishable from the Airy distribution of Eq. (7c).
In the literature, transmittance responses are commonly classified based on the finesse parameter, whose value $4\mathcal {R}/(1-\mathcal {R})^2$ increases with the reflectivity. For low finesse devices, the response closely resembles a pure sinusoid. This characteristic allows for increased throughput, resulting in higher signal to noise ratio (SNR) captured by the sensors. In the case of high finesse, the spectral response exhibits sharper peaks, resulting in an enhanced periodic bandpass filtering effect. Different finesse regimes can be determined for each wave model, by establishing the maximum reflectivity value such that a given error measure between the transmittance response of the $W$-wave model and the Airy distribution remains below a certain threshold. Figure 5(b) illustrates this concept, employing the root mean square error (RMSE) as the chosen error measure.
2.3 Proposed formulation of the image formation model
In order to characterize the overall spectral response of the instrument at a given pixel, the physical acquisition model employed from Eq. (2) may be simplified, assuming that the optical transmittance response is roughly constant within the targeted solid angle $\Omega _j$ and $\cos \theta \approx 1$:
Here, $T_{\boldsymbol{\beta }}(\sigma )$ models the transmittance response of the entire instrument associated to a given pixel on the focal plane array (FPA). For convenience, it is useful to describe it in terms of the expression of Eq. (8):
The terms $\mathcal {A}(\sigma )$ and $\mathcal {R}(\sigma )$ exhibit strong coupling in Eq. (11). To estimate their contributions separately, we impose them to be slowly varying functions with respect to $\sigma$. To achieve this, we restrict their models to polynomials of limited degree $N_d$:
The OPD value $\delta$ is assumed to be constant with the wavelength $\sigma$ as the rays interfer within the air in the prototypes under test following the design of Fig. 1(a). This assumption is extended to the phase shift $\varphi ^{}_0$ in order to simplify the computation. These last hypotheses may be too limiting for interferometric cavities made of dispersive materials. For such cases, one may suppose a prior knowledge of this dispersion as function of $\sigma$ to reduce the problem to the estimation of the interferometer thickness, which is independent on $\sigma$.
Our goal then summarizes to find an estimation $\hat {\boldsymbol{\beta }}$ of the $2N_d+4$ elements of $\boldsymbol{\beta }=\left [a_0\,,\ldots,\,a_{N_d},\,r_0\,,\ldots,\,r_{N_d},\,\delta,\,\varphi _0^{}\right ]$ which allows to approximate the transmittance response $T_{\boldsymbol{\beta }}$ as accurately as possible.
3. Proposed characterization procedure
In this section, we present the proposed procedure for the spectral characterization of FP interferometers. Specifically, we describe the measurement setup for the characterization in Section 3.1 and we provide an overview of IRCA in Section 3.2, detailing each of its composing steps in the subsequent sections.
3.1 Measurement setup for the characterization of the device
To accurately characterize a given device under test and allow the inference of its parameters, it is necessary to capture a specific set of observations from reference sources under controlled conditions. Perhaps the most straightforward approach involves illuminating the device with a flat field illumination using a set of monochromatic incident spectra with predefined central wavenumbers. In fact, based on Eq. (9), the corresponding $N_a$ acquisitions $\mathbf {y}\in \mathbb {R}^{N_a}$ are samples of the expected value of $T_{\boldsymbol{\beta }}(\sigma )$ evaluated at the specific wavenumbers $\boldsymbol{\sigma }\in \mathbb {R}^{N_a}$.
In order to achieve this result, we propose the measurement setup shown in Fig. 6. It involves the utilization of a wideband lamp as the light source, whose emitted light is filtered by a monochromator (i.e., equipped with a diffraction grating). The bandwidth of the monochromator is deliberately narrower than the spectral resolution of the device, ensuring a sharply impulsive filtered spectrum. Subsequently, this spectrum is uniformly scattered over the device under test by means of an integrating sphere. By tuning the monochromator, a series of central wavenumbers is selected, and the device under test captures an image for each illumination in sequence.
Fig. 6. Measurement setup for the characterization of multi-aperture interferometric imaging spectrometers.
An external spectrometer or probe is used to measure the incident power of the instrument. The measured value is used to equalize the energy of all the acquired images at different wavenumbers, with background level set to zero. Finally, the vector $\mathbf {y}$ is obtained by extracting the specific spatial position from the acquired datacube, corresponding to the pixel being characterized.
Therefore, we formalize the problem at hand as finding the estimation $\hat {\boldsymbol{\beta }}$ of the parameter vector such that:
3.2 Overview of the interferometer response characterization algorithm (IRCA)
Solving Eq. (13) is a particularly challenging problem, due to the nonlinear dependency of $T_{\boldsymbol{\beta }}$ from the parameters $\boldsymbol{\beta }$. The available tools for solving nonlinear regression methods are particularly sensitive to converging to non-local maxima [24], so that a proper initialization is critical to produce an accurate parametrization of the optical system.
The proposed interferometer response characterization algorithm (IRCA), depicted in Fig. 7, consists of three steps, each dedicated to processing one of the three different sufficient statistics extracted from the $N_a$ images captured during the measurement session. This approach is designed to enhance the overall robustness of the final result by leveraging multiple aspects of the available information. We describe each of these steps below:
- • The gain estimation step processes the vector $\mathbf {w}\in {\mathbb {R}^{N_a}}$, which represents the flat field statistic. This vector is used to obtain an initial assessment of the gain coefficients $\{\hat {a}_i\}_{ {i} \in [{0}, \ldots , {N_d}]}$ of $\mathcal {A}(\sigma )$. The flat field statistic captures the response of the pixel under test without the presence of the interferometric fringes. In cases where the vector $\mathbf {w}$ is not directly available, it can be approximated by evaluating the percentile from the raw acquisition across the entire focal plane. This approach takes advantage of the global response of the image, which naturally dampens the oscillations caused by the interferometric fringes.
- • The maximum likelihood (ML) initialization involves processing the vector $\mathbf {u}\in \mathbb {R}^{N_a}$, which represents the mean over neighbors. This step returns an initial estimation $\left [\hat {\delta },\,\hat {r}_0,\,\hat {\varphi }_0\right ]$ of the remaining parameters, namely the OPD, reflectivity, and phase shift, respectively. The vector $\mathbf {u}$ is obtained by calculating the spatial average of the raw acquisition within a square window centered around the pixel under test. This averaging helps reduce the noise associated with the acquisition and could be performed even in the temporal domain, if such information is available. At this stage of the estimation process, the parameters are assumed to be constant with the wavenumbers.
- • In the trust region refinement (TRR) step, we process the raw acquisition $\mathbf {y}\in \mathbb {R}^{N_a}$ to obtain the final estimation $\hat {\boldsymbol{\beta }}$ of the complete set of parameters. To achieve this, we initialize a LM algorithm [21] with the parameter vector $\boldsymbol{\beta }^{[0]}$ whose elements were inferred in the previous steps. Subsequently, we iterate through the algorithm to solve Eq. (13).
Fig. 7. Overview of the proposed IRCA algorithm for the characterization of the transmittance response of a single Fabry-Perot interferometer.
The following sections describe each of these steps in further detail. It is important to note that in certain scenarios, such as non-imaging systems where a one-dimensional acquisition is obtained using a single pixel sensor, the mean over neighbors and flat field statistic may not be available. Nonetheless, in these cases, the algorithm can still be applied by setting $\mathbf {u}$ equal to $\mathbf {y}$ and the elements of $\mathbf {w}$ equal to the average value of $\mathbf {y}$.
3.3 Step 1: gain estimation
We formalize the problem associated to the gain estimation as follows:
where $\hat {\mathbf {a}}=\{\hat {a}_m\}_{ {m} \in [{0}, \ldots , {N_d}]}$ describes our estimation the coefficients of the polynomial representation $\hat {\mathcal {A}}(\sigma )=\sum _{m=0}^{N_d}\hat {a}_m\sigma ^m$ of the gain $\mathcal {A}(\sigma )$. The vector $\mathbf {w}=\{w_i\}_{ {i} \in [{1}, \ldots , {N_a}]}$ contains the samples of the flat field statistic.We propose to solve the problem above with a nonlinear regression approach using the LM algorithm with the implementation of [21], as described in Appendix A. We initialize the gain coefficients vector $\mathbf {a}$ by setting the first element to the average value of $\mathbf {w}$ and the rest to zero.
3.4 Step 2: maximum likelihood (ML) initialization
The maximum likelihood (ML) initialization step defines a procedure that is as an extension of the simplistic OPD estimation algorithm proposed in [17]. In the step, we assume to operate in a low finesse regime, so that the transmittance response $T_{\boldsymbol{\beta }}(\sigma )$ behaves like the 2 waves model of Eq. (7a):
In the above equation, we implicitly impose that the reflectivity $\mathcal {R}$ is uniform and equal to $r_0$ over the whole wavenumber range. By normalizing both terms of the minimization problem of Eq. (13) and applying it to the mean over neighborhood vector $\mathbf {u}$, the problem can be rewritten as:
Equation (16) is in the form of the classical problem of the inference of the parameters in a sinusoid affected by Gaussian noise, which is a well known problem in the literature [25, Example 7.16]. Specifically, it is a well known result that the maximum likelihood estimator $\hat {\delta }$ is equal to the OPD value which maximizes the periodogram:
In other terms, the estimator is the value that maximizes the generalized discrete Fourier transform (DFT) of $\mathbf {v}$. The above result is valid for values of $\delta$ reasonably far from the extremes of the interval $[0, 1/(2\Delta \sigma )]$, where $\Delta \sigma =(\sigma _{max}-\sigma _{min})/N_a$ is the average central wavenumber step.
Equation (17) can be solved numerically over a sampled version of the interval of interest, yet the accuracy is limited to a resolution of $1/(2N_a\Delta \sigma )$. If the OPD is approximately known, it is computationally efficient to evaluate Eq. (17) within a reduced interval centered around its nominal value.
The estimation $\hat {r}_0$ of the reflectivity $r_0$ is then obtained in terms of the the estimation $\hat {\alpha }$ of the amplitude $\alpha$ of the sinusoid:
In the above equation, $\arctan$ denotes the four-quadrant arctangent version that allows for $\hat {\varphi }_0^{}$ to assume any value in the range $[-\pi,\pi )$. The ML method requires very low computational power, but its applicability is limited by the validity of its assumptions. Some other possible initialization strategies, such as the exhaustive search (ES) developed in [26] which is based on a grid search in the sample space of the parameters, have the advantage to work with a wider variety of models. They are however vastly slower and may not necessarily produce more accurate results, as the estimations for $r_0^{}$ and $\varphi _0^{}$ are limited to the finite amount of values of the discrete sample space.
3.5 Step 3: trust region refinement (TRR)
The final parameter estimation $\hat {\boldsymbol{\beta }}$ follows a similar procedure as described in Section 3.3. Specifically, the LM algorithm is employed once again, but this time to solve Eq. (13). The parameter vector $\boldsymbol{\beta }$ is initialized with values from $\boldsymbol{\beta }^{[0]}=[\hat {a}_0,\ldots,\,\hat {a}_{N_d},\, \hat {r}_0^{},\, 0,\ldots,\,0,\, \hat {\delta },\hat {\varphi }_0^{}]$, where the non-zero elements correspond to the coefficients estimated at the step 1 and 2 of the algorithm.
4. Experimental results
This section presents the experimental results obtained from the characterization of a series of ImSPOC prototypes with various characteristics. In Section 4.1 we describe the experimental setup, in Section 4.2 we test various configurations for the proposed algorithm and compare its performances with previous works. Finally, in Section 4.3 we discuss the physical interpretation of the parameters. A Python implementation of the proposed algorithms, together with a simulator of the image formation for multi-aperture Fabry-Perot imaging spectrometers, is available at the first author’s repository [27].
4.1 Experimental setup
For this work, the characterization datacubes were captured with the setup shown in Fig. 6, using a tunable monochromatic light source from Zolix Instruments Co., Ltd, with a 500 W Xenon light source model Gloria-X500A and a monochromator model Omni-300$\lambda$i. We also utilized a 5.3-inch diameter integrating sphere coated with Spectralon (model 4P-GPS-053-SF from Labsphere, Inc.). The incident optical power was measured either with the fiber optic gated spectrometers model USB2000+ from Ocean Optics, Inc. or with the photodiode power sensor model S120VC from Thorlabs. The product specifications can be found on the websites of the respective manufacturers.
The devices under test are four different ImSPOC prototypes, whose characteristics are described in Table 2. Each prototype features an array of interferometers disposed over a bidimensional matrix in a staircase pattern, whose thicknesses linearly increase with a nominally constant step size $\Delta d$. While sharing the same underlying concept, each prototype is specifically designed for different applications. Prototype 1 and 2 are specifically tailored for the measure of atmospheric pollution, prototype 3 functions as an imaging system for capturing the phenomenon of northern lights, and prototype 4 is intended for greenhouse gas detection. In terms of spectral sensitivity, prototypes 1 to 3 operate within the visible/ultraviolet wavelength range, whereas prototype 4 covers the near-infrared spectrum. For each device, a characterization datacube was captured using the procedure described in Section 3.1. The central wavenumbers of the monochromator are chosen to be regularly spaced with a step size $\Delta \sigma$. The wavenumber step size $\Delta \sigma$ is selected to satisfy the Nyquist condition $\Delta \sigma < 1/(2\delta_{max})$. This selection ensures that aliasing effects are avoided in the sampling of the transmittance response of the interferometer with the largest OPD $\delta_{\max }=2\left(N_i-1\right) \Delta d$. The condition is satisfied in all experimental setups, although only by a small margin for prototype 1 and to a lesser extent for prototype 3, due to time constraints. The specifications for these measurements are reported in Table 2.
Table 2. Characteristics of the available ImSPOC prototypes used in this work and of their spectral characterization experimental acquisitions.
Given a characterization datacube and a specific subimage within it, the central pixel of the chosen subimage is extracted to construct the raw acquisition vector $\mathbf {y}$. The mean over neighbors is computed using a $11 \times 11$ kernel window centered around the extracted pixel, while the 90-percentile metric is instead employed as the flat field statistic. Next, we apply the characterization method described in Section 3 to obtain the characterization vector $\hat {\boldsymbol{\beta }}$, with $N_d=5$ as the degree of the polynomial for the reflectivity and gain. To verify the quality of the estimation, we use the RMSE metric, defined as follows:
4.2 Algorithm and model comparisons
The IRCA is tested here with different configurations, and we assume that the gain estimation is always carried as a pre-processing step. We compare its results with previous works [17,26] which we can conveniently frame within our proposed framework.
We employ different wave models for the optical transmittance response $T_{\boldsymbol{\beta }}(\sigma )$, according to the definitions of Eqs. (7) and (10), specifically for the case of $W=2$, $3$, or $\infty$ emerging light rays.
The RMSE results, given in Table 3, shows that, when all the three steps are performed, the proposed method is consistently the best performing, regardless of the different characteristics of the prototypes. It also highlights how the $\infty$-wave model, which is a better representation of the generalized Airy distribution, provides a more accurate fit for the spectral response. Both tested initializations reach comparable results, suggesting that the ML method should be preferred, as it is faster by a factor of 10-20 times over the ES methodology.
Table 3. Model characterization RMSE comparison. Best results are in bold.
The proposed procedure was tested both with and without the trust region refinement step, in order to showcase the advantage of the iterative curve fitting procedure. The TRR step has a considerable impact on the accuracy of the results. This is due to the LM algorithm exploring a continuous space of parameters, resulting in better performance with respect to the ML by itself where the OPD space is explored in discrete steps. A visual comparison between their reconstructed spectral responses is shown in Fig. 8. While the ML algorithm by itself already infers the OPD with a remarkable accuracy, the curve does not follow accurately the trend of the data, as the other parameters vary significantly with the wavenumbers.
Fig. 8. Spectral responses of the characterization results with the IRCA. (a), (b) Two examples of fitting the acquisition with the transmittance response using estimated parameters. The acquisition (in orange) is compared to the estimation with the proposed IRCA method (in green) and with the ML method (in blue). (c) Estimated gain. In the orange curve, we allow the LM algorithm to freely update the parameters of the gain polynomial, while in the blue curve we only allow to change the scaling factor of the initialized gain. (d) Estimated reflectivity. The shaded area around the curves represents one standard deviation interval across different interferometers.
The results for the prototype 2 for the estimation of the gain and reflectivity are given in Fig. 8. The figures showcase the advantage of the TRR step in assessing the gain for each detector separately. The analysis of the reflectivity reveals a reduced sensitivity of the instrument at the extreme values of the wavenumber range, which aligns with the expected spectral response of the reflective coating.
4.3 Physical interpretation of the results
For the prototype 2, the measurements for the characterization was repeated at three different dates, using progressively refined setups. The proposed method was applied to each of those datasets in order both to analyze the robustness of the algorithm and to detect eventual drifting in the parameters. Figure 9(a) provides a visual comparison of the results. For the normal incidence of the light illumination, the OPDs are roughly expected to be double the thickness of interferometers, as a consequence of substituting $\theta =0$ in Eq. (5). The estimation of the OPDs in Fig. 9(a) stays vastly consistent across all sessions, with only very sparse examples where the results do not align.
Fig. 9. Characterization of prototype 2 from measurements taken at different times. (a) Visualization of the estimated OPDs for all the interferometers, arranged in ascending order of nominal thicknesses. The yellow halo indicates the region that was not explored for the ML initialization, as that strays too far from the nominal values known from the design of the instrument. (b), (c) Comparison between the estimated transmittance responses (blue) and raw acquired data (orange) for interferometer #150. There are significant differences in the intensity values and appearance of the orange curves between the two figures. This stems from the fact that the raw acquisitions on 2021-12-23 were divided by the spectral response of the light source, while the acquisitions on 2021-10-20 are unaltered as we lack this information for that particular session.
The figure shows a pattern of alternating slopes. To make sense of this effect, we also present the estimation of the OPDs for the prototype 1. Prototype 1 is designed with a distinctive staircase pattern, characterized by cavity thicknesses that increase in constant steps on both sides of a central vertical axis. This pattern is reflected in the deviation of the estimated OPD increase with respect to the expected value. As depicted in Fig. 10, the observed discrepancy can likely be attributed to a tilt of the optical plates relative to the desired parallel installation.
Fig. 10. Estimated OPDs for the prototype 1. (a) The heatmap presents the variation in the increase of the OPD between successive interferometers. The indexes of the interferometers are arranged in ascending order of their nominal thickness. (b), (c) Effect of tilting the optical plate on with respect to an aligned one. The tilt can result in either a compression (indicated by red shades) or an expansion (indicated by blue shades) of the air gap in the cavity with respect to the nominal value. The indices inside the cavities denote the positions of the interferometers within the row.
We also employed the proposed algorithm for the spectral characterization of prototype 2 at different angles of incidence, following the common practice in the literature [28–30]. Specifically, we conducted a comprehensive characterization of the transmittance response for every pixel on the focal plane using the ML+TRR variant of the proposed IRCA method. The resulting estimation for the OPD is displayed in Fig. 11, where the subimage area of the instrument corresponds to a field of view of around $\pm 5$ degrees. Our measurements indicate that the relative variation of the OPD with respect to its value in the central pixel is experimentally found to be more or less identical across different interferometers. The analysis shows a decreasing trend for the OPD which emanates radially from a central position. This is an expected result from Eq. (5), as the OPD decreases when the angle of incidence increases. The true optical axis is however shifted to the left with respect to the geometrical center that we have arbitrarily chosen as center of the subimage, as shown in Fig. 11(b). On the contrary, the estimated values of the reflectivity preserve a certain flatness within their own subimage, as depicted in Fig. 11(c). The spatial analysis also allows for the detection of certain instrument defects, such as those of two of the subimages in bottom left area of Fig. 11(a). In a significant portion of these subimages, the device exhibits behavior that significantly deviates from the expected one, suggesting the presence of defects in the reflective coating.
Fig. 11. Characterization of the prototype 2 over the entire field of view of the device using the IRCA method. (a) Illustration of the relative increase of the estimated OPD compared to its value at the center of the corresponding subimage, indicated by a cross. Red and blue shades indicate values larger and smaller than the reference, respectively. Black pixels correspond to areas where the acquisition data was unavailable or where the LM algorithm did not converge within 100 iterations. (b) Zoomed-in visualization of the yellow framed area. (c) Estimated reflectivity per pixel over the same zoomed area, averaged over the wavenumber range of the device.
5. Conclusion
In this paper, we presented a characterization procedure for multi-aperture imaging spectrometers based on Fabry-Perot interferometers. We described the image formation model and we expressed its transmittance response in terms of a limited amount of parameters, following the formulation of Airy distribution and describing different regimes of finesse under the same framework. The proposed characterization algorithm exploits the two emerging wave model to cast the problem in the Fourier domain, where the ML estimator for the OPDs is equivalent to a maximization of the periodogram. This result is then refined through nonlinear regression. Using the proposed multiple-step procedure in the algorithm allows for both robustness and precision in the final results. The estimated parameters can highlight manufacturing issues in an easily interpretable format. A proper characterization is extremely important for a proper recovery of the spectrum, which in the future could be optimized jointly with the spectral response of the devices in architectures where such parameters can be learned dynamically [31,32].
A. Levenberg-Marqardt algorithm
We aim to provide here an approachable explanation of the LM algorithm, to build an intuition of what are the operations involved in the process. The algorithm can be seen as a trust region based approach for nonlinear regression. The aim is to find an estimation $\hat {\boldsymbol{\beta }}$ of the parameters $\boldsymbol{\beta }=\{\beta _m\}_{ {i} \in [{1}, \ldots , {N_m}]}$, in order for the samples $\{t_i(\boldsymbol{\beta })\}_{ {i} \in [{1}, \ldots , {N_a}]}$ of an analytical function to fit a set of observation $\{y_i\}_{ {i} \in [{1}, \ldots , {N_a}]}$.
The algorithm addresses the problem by finding a sequence of iteratively more accurate solutions $\{\boldsymbol{\beta }^{[q]}\}_{q\ge 0}$ from a given initialization $\boldsymbol{\beta }^{[0]}$, using the following update rule:
Equation (21) admits as analytical solution:
Funding
AuRA region and FEDER (ImSPOC-UV: convention FEDER n. RA0022348); Agence Nationale de la Recherche (FuMultiSPOC: ANR-20-ASTR-0006, LabEx FOCUS: ANR-11-LABX-0013).
Acknowledgments
The authors would like to thank the anonymous reviewers for their constructive feedback.
For illustration, authors Yann Ferrec, and Etienne Le Coarer, are represented below as YF, and ELC
Disclosures
YF, ELC: ONERA, UGA - U.S. Patent: 10,677,650 B2 (P).
Data Availability
Data underlying the results presented in this paper are available in Ref. [27].
References
1. A. K. Tilling, G. J. O’Leary, J. G. Ferwerda, S. D. Jones, G. J. Fitzgerald, D. Rodriguez, and R. Belford, “Remote sensing of nitrogen and water stress in wheat,” Field Crop. Res. 104(1-3), 77–85 (2007). [CrossRef]
2. E. Ben-Dor, S. Chabrillat, J. A. M. Demattê, G. R. Taylor, J. Hill, M. L. Whiting, and S. Sommer, “Using imaging spectroscopy to study soil properties,” Remote. Sens. Environ. 113, S38–S55 (2009). [CrossRef]
3. E. Adam, O. Mutanga, and D. Rugege, “Multispectral and hyperspectral remote sensing for identification and mapping of wetland vegetation: A review,” Wetl. Ecol. Manag. 18(3), 281–296 (2010). [CrossRef]
4. Y. C. Kim, H. G. Yu, J. H. Lee, D. J. Park, and H. W. Nam, “Hazardous gas detection for FTIR-based hyperspectral imaging system using DNN and CNN,” Proc. SPIE 10433, 1043317 (2017). [CrossRef]
5. S. Gousset, L. Croizé, E. Le Coarer, Y. Ferrec, J. Rodrigo Rodrigo, and L. Brooker, “NanoCarb hyperspectral sensor: On performance optimization and analysis for greenhouse gas monitoring from a constellation of small satellites,” CEAS Space J. 11(4), 507–524 (2019). [CrossRef]
6. R. A. Borsoi, T. Imbiriba, J. C. M. Bermudez, C. Richard, J. Chanussot, L. Drumetz, J.-Y. Tourneret, A. Zare, and C. Jutten, “Spectral variability in hyperspectral data unmixing: A comprehensive review,” IEEE Geosci. Remote Sens. Mag. 9(4), 223–270 (2021). [CrossRef]
7. J. Bacca, E. Martinez, and H. Arguello, “Computational spectral imaging: a contemporary overview,” J. Opt. Soc. Am. A 40(4), C115–C125 (2023). [CrossRef]
8. L. Huang, R. Luo, X. Liu, and X. Hao, “Spectral imaging with deep learning,” Light: Sci. Appl. 11(1), 61 (2022). [CrossRef]
9. A. Dolet, D. Picone, S. Gousset, M. Dalla Mura, E. Le Coarer, and D. Voisin, “Using zenith observations for evaluation of an improved interferometric imaging spectrometer,” in Proceedings of European Geosciences Union General Assembly (EGU, 2021), paper EGU21-2536.
10. Y. Oiknine, I. August, and A. Stern, “Multi-aperture snapshot compressive hyperspectral camera,” Opt. Lett. 43(20), 5042–5045 (2018). [CrossRef]
11. Y. Ferrec, G. Bonnery, L. Brooker, L. Croizé, S. Gousset, and E. Le Coarer, “NanoCarb part 1: Compact snapshot imaging interferometer for CO2 monitoring from space,” Proc. SPIE 11180, 1118021 (2019). [CrossRef]
12. S. Gousset, L. Croizé, E. Le Coarer, Y. Ferrec, and L. Brooker, “NanoCarb part 2: Performance assessment for total column CO2 monitoring from a nano-satellite,” Proc. SPIE 11180, 111803Q (2019). [CrossRef]
13. S. Gousset, N. Guerineau, L. Croizé, E. Le Coarer, T. Laveille, and Y. Ferrec, “NANOCARB-21: a miniature Fourier-transform spectro-imaging concept for a daily monitoring of greenhouse gas concentration on the Earth surface,” Proc. SPIE 10562, 105624U (2017). [CrossRef]
14. N. Guerineau, E. Le Coarer, Y. Ferrec, and F. De La Barriere, Fourier transform multi-channel spectral imager, U.S. Patent: 10,677,650 B2 (9 June 2020).
15. M. Jouni, D. Picone, and M. Dalla Mura, “Model-based spectral reconstruction of interferometric acquisitions,” in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (IEEE, 2023), (to be published).
16. D. Picone, “Model based signal processing techniques for nonconventional optical imaging systems,” Ph.D. dissertation (Université Grenoble Alpes, 2021).
17. A. Dolet, D. Picone, M. Dalla Mura, D. Voisin, S. Gousset, S. Douté, and E. Le Coarer, “Gas characterization based on a snapshot interferometric imaging spectrometer,” Proc. SPIE 11155, 1115502 (2019). [CrossRef]
18. Y. Oiknine, I. August, D. G. Blumberg, and A. Stern, “NIR hyperspectral compressive imager based on a modified Fabry-Perot resonator,” J. Opt. 20(4), 044011 (2018). [CrossRef]
19. M. Pisani and M. E. Zucco, “Compact imaging spectrometer combining Fourier transform spectroscopy with a Fabry-Perot interferometer,” Opt. Express 17(10), 8319–8331 (2009). [CrossRef]
20. M. Zucco, M. Pisani, V. Caricato, and A. Egidi, “A hyperspectral imager based on a Fabry-Perot interferometer with dielectric mirrors,” Opt. Express 22(2), 1824–1834 (2014). [CrossRef]
21. J. J. Moré, “The Levenberg-Marquardt algorithm: implementation and theory,” in Numerical analysis, (Springer, 1978), pp. 105–116.
22. U. C. Hasar, I. Y. Ozbek, and T. Karacali, “Nondestructive optical characterization of Fabry-Pérot cavities by full spectra fitting method,” IEEE Photonics Technol. Lett. 30(15), 1404–1407 (2018). [CrossRef]
23. N. Ismail, C. C. Kores, D. Geskus, and M. Pollnau, “Fabry-Pérot resonator: spectral line shapes, generic and related Airy distributions, linewidths, finesses, and performance at low or frequency-dependent reflectivity,” Opt. Express 24(15), 16366–16389 (2016). [CrossRef]
24. A. P. Ruszczyński, Nonlinear optimization (Princeton University, 2006).
25. S. M. Kay, Fundamentals of statistical processing (Prentice Hall, 1993). Vol. 1, pp. 193–195.
26. D. Picone, A. Dolet, S. Gousset, D. Voisin, M. Dalla Mura, and E. Le Coarer, “Characterisation of a snapshot Fourier transform imaging spectrometer based on an array of Fabry-Perot interferometers,” in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (IEEE, 2020), pp. 1529–1533.
27. D. Picone, “Code repository for the interferometer response characterization algorithm (IRCA),” Zenodo (2023) [Retrieved 1 June 2023], http://dx.doi.org/10.5281/zenodo.7978514.
28. W. Zhang, H. Song, X. He, L. Huang, X. Zhang, J. Zheng, W. Shen, X. Hao, and X. Liu, “Deeply learned broadband encoding stochastic hyperspectral imaging,” Light: Sci. Appl. 10(1), 108 (2021). [CrossRef]
29. S. Feng, Z. Wang, X. Cheng, and X. Dun, “Superposition Fabry-Perot filter array for a computational hyperspectral camera,” Opt. Lett. 48(5), 1156–1159 (2023). [CrossRef]
30. J. Yang, K. Cui, Y. Huang, W. Zhang, X. Feng, and F. Liu, “Angle-insensitive spectral imaging based on topology-optimized plasmonic metasurfaces,” arXiv, arXiv:2212.07813 (2022). [CrossRef]
31. Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” IEEE Trans. Signal Process. 65(3), 794–816 (2016). [CrossRef]
32. Y. Yang, J. Sun, H. Li, and Z. Xu, “ADMM-CSNet: A deep learning approach for image compressive sensing,” IEEE Trans. Pattern Anal. Mach. Intell. 42(3), 521–538 (2020). [CrossRef]
