Fourier multispectral imaging

Current multispectral imaging systems use narrowband filters to capture the spectral content of a scene, which necessitates different filters to be designed for each application. In this paper, we demonstrate the concept of Fourier multispectral imaging which uses filters with sinusoidally varying transmittance. We designed and built these filters employing a single-cavity resonance, and made spectral measurements with a multispectral LED array. The measurements show that spectral features such as transmission and absorption peaks are preserved with this technique, which makes it a versatile technique than narrowband filters for a wide range of multispectral imaging applications. © 2015 Optical Society of America OCIS codes: (110.4234) Multispectral and hyperspectral imaging; (240.0310) Thin films. References and links 1. M. T. Eismann, Hyperspectral remote sensing (Bellingham: SPIE Press, 2012). 2. B. Geelen, N. Gack, and A. Lambrechts, “A compact snapshot multispectral imager with a monolithically integrated per-pixel filter mosaic,” in SPIE MOEMS-MEMS, International Society for Optics and Photonics, 89740L89740L (2014). 3. D. Yi and L. Kong, “Fabrication of densely patterned micro-arrayed multichannel optical filter mosaic,” Journal of Micro/Nanolithography, MEMS, and MOEMS 10(3), 033020-033020 (2011). 4. G. Themelis, J. S. Yoo and V. Ntziachristos, “Multispectral imaging using multiple-bandpass filters,” Opt. Lett. 33, 1023-1025 (2008). 5. K. C. Balram and D. A. B. Miller, “Self-aligned silicon fins in metallic slits as a platform for planar wavelengthselective nanoscale resonant photodetectors,” Opt. Express 20, 22735-22742 (2012). 6. K. Hirakawa and K. J. Barnard, “Fourier spectral filter array design for multispectral image recovery,” in OSA Imaging Systems and Applications, Optical Society of America (2014), pp. IM1C-5. 7. J. Jia and K. Hirakawa, “Single-shot fourier yransform multispectroscopy,” in IEEE International Conference on Image Processing, Institute of Electrical and Electronics Engineering (2015). 8. C. Ni, J. Jia, K. Hirakawa, A. M. Sarangan, “Design and fabrication of sinusoidal spectral filters for multispectral imaging,” in SPIE Optics + Photonics, International Society for Optics and Photonics (2015). 9. J. W. Goodman and R. L. Haupt, Statistical optics (John Wiley & Sons, 2015). 10. A. Chakrabarti and T. Zickler, “Statistics of real-world hyperspectral images,” in Proc. IEEE Conf. on Computer Vision and Pattern Recognition (2011), pp. 193-200. 11. H. A. Macleod, Thin-film optical filters (Institute of Physics Publishing Bristol and Philadelphia, 2001).


Introduction
Conventional colorimetric image sensors summarize the spectral content of the image by only three values metameric to human visual system.By contrast, hyperspectral imaging (HSI) makes dense measurements over a specified wavelength range.The HSI data is useful for making complete spectral characterization of objects, providing the positional information of objects that can be used for the higher level process such as tracking and object recognition.A third alternative is multispectral imaging (MSI).Instead of sampling densely over a continuous spectra, a finite number of spectral measurements are made by N predetermined narrowband filters.This significantly reduces the hardware complexity, albeit at the loss of some spectral density.However, we demonstrate in this paper that the spectral features can be recovered to a great extent by using unconventional filtering techniques instead of simple band-pass elements.
Existing multispectral imaging modalities capture spectral content of a scene using narrowband filters (or their combinations) that attenuate light outside specific wavelength ranges [1][2][3][4].The center wavelengths of the pass bands are designed to match the object's spectral features of interest such as emission, reflection, or absorption peaks that arise from their material and chemical composition [5].However, these are application-specific-different spectral features and applications would require different filters to be synthesized.Hence these narrowbandbased MSI modalities do not generalize well to multi-purpose MSI, where multiple spectral peaks and absorption bands need to be observed, or when we lack a priori knowledge of the material composition.
We addressed this shortcoming by developing an alternative to the narrowband filter-based MSI modality.Our goal is to find a better-suited basis set for spectral decomposition by choosing filter transmission functions that are fundamentally different than the bandpass filter.Specifically, we demonstrate an innovative technique which utilizes multichroic filters that have sinusoidal transmission spectra in the wavenumber domain.Unlike narrowband measurements, multichroic filter measurements largely avoid aliasing that contaminates spectra when undersampling with narrowband filters.Because of this, there are a number of advantages to multichroic filters with different sinusoidal transmission responses, including the better spectrum recovery from discrete filters, the preservation of spectral features such as peaks and absorption bands that are critical for target detection, and the direct correspondence to single-cavity resonator implementation.Due to its foundation in Fourier transform spectroscopy, we refer to the proposed multichroic filter-based image acquisition technique as Fourier MSI.
The use of sinusoidal transmission filters was first proposed in [6][7][8].However, previous reporting lacked the experimental validations and direct correspondences to multichroic filters that are the focus of this work.We employed dielectric interference filters that afford a greater performance and spectral flexibility compared to polymeric dye-doped absorptive filters.Multilayer interference filters are based on reflections at the dielectric interfaces and phase shifts within the layers.By carefully choosing the type, sequence, and thicknesses of the thin film materials, any desired spectral profile can be achieved, while reflecting all out-of-band light.

Fourier spectral filter
Let us first review the Fourier transform spectroscopy.Typical spectral imagery is represented as a cube with two dimensions in space and one spectral dimension in wavelength or wavenumber.Assuming wavenumber is the unit of choice for the spectral dimension, let X : Z 2 ×R → R + be an image cube, where X(n n n, σ ) denotes the spectral radiance measured at the pixel location n n n = (n 0 , n 1 ) T ∈ Z 2 and wavenumber σ ∈ R. Invoking the Wiener-Khinchin theorem, performing an inverse Fourier transform on X(n n n, σ ) with respect to σ yields an autocorrelation Intuitively, OPD ζ refers to a coherent length, typically expressed in the units of nm (since it is the reciprocal to wavenumber σ in nm −1 ).Fig. 1 illustrates an example hyperspectral cube in spatial-OPD (X(n n n, ζ )) domain.
In this paper, we propose undersampling X(n n n, ζ ) in the spatial-OPD domain as an attractive general purpose MSI modality.That is to say, we propose to represent image cube X(n n n, σ ) with 0nm 420nm 840nm 1260nm 1680nm 2100nm Fig. 1.Hyperspectral image cube in Space-OPD representation.Shown using dataset in [10].
only K OPD samples {Y 0 (n n n), . . .,Y K−1 (n n n)} where for some sampling interval ∆ζ > 0. We refer to the OPD samples {Y 0 (n n n), . . .,Y K−1 (n n n)} collectively as Fourier multispectral imaging, or Fourier MSI.As we prove below, Fourier MSI has clear advantages over the conventional narrowband filtering approaches to MSI.
Recall that the quantum efficiency of a detector is limited to a spectral range (e.g.σ ∈ (1000nm −1 , 380nm −1 ) in CMOS sensors).It implies that discrete OPD samples X(n n n, k/2σ MAX ) for integer k ∈ Z suffices for reconstructing X(n n n, σ ) in the ranges σ ∈ (0, σ MAX ): where ∆ξ = 1 2σ MAX .Even fewer samples are needed to represent the spectra X if the detector only responds to the range σ ∈ ( 1 2 σ MAX , σ MAX )-in order for Eq. ( 3) to be valid, ∆ξ = 1 σ MAX is sufficient.
We conclude that the dense or continuous sampling of X(n n n, ζ ) in the OPD domain is unnecessary.The OPD sampling interval of ∆ζ = 1 σ MAX in Eq. ( 3) suffices for many applications.Fig. 1 illustrates hyperspectral image signal in Space-OPD representation.
To further reduce the number of spatial-OPD measurements X(n n n, ζ ), consider keeping only the first K OPD samples, as Eq. ( 2) suggests.There are three main reasons for why , meaning it is unnecessary to acquire negative OPD samples.Second, OPD representation enjoys energy compaction.This is due in part to the fact that ζ = 0 is the global maximum for any autocorrelation function: This observation is also consistent with the well-accepted fact that the coherent length of a typical waveform decays rapidly.Fig. 1 further confirms this trend.
Third, recall that emission, reflection, or absorption peaks are key spectral features for material/chemical identification applications.Fourier MSI in Eq. ( 2) guarantees that the locations of these spectral features are well preserved.To see why this is the case, consider spectrum where I K (k) is an indicator function that activates when |k| < K, and I K (σ ) is its discrete time Fourier transform.Noting that I K (σ ) is a sinc function, the convolution in Eq. ( 5) has the effect of broadening the spectral emission/reflection/absorptions peaks.Thanks to the symmetry of sinc function, the location of the spectral peaks in X(n n n, σ ) are preserved, even when K is small.The advantage to increasing the number of filters (larger K) is that the peaks close to each other can be distinguished (because the widths of the sinc function decreases).
Based on the analysis above, we are interested to make OPD measurements {Y 0 (n n n), . . .,Y K−1 (n n n)}.Known as Michaelson-Morley interferometry, one such way is to combines two optical fields at the detector, one delayed in time/optical path with respect to the other.However, this classical Fourier transform spectroscopy setup does not preserve pixel structures.We propose to make OPD measurements using optical spectral filter with specific sinusoidal spectral transmission functions, resulting in multispectral measurements of the form where α k is the offset and χ k is the so-called filter "contrast"-a notion we define more rigorously in Eq. ( 9) below.The desired OPD samples can be recovered from the filtered measurements:

Design of Fourier spectral filters
We now develop filters of the Eq. ( 6).Most conventional dielectric filters are designed for sharp spectral transitions with high contrast, like high-pass, low-pass, notch, etc. whose features emerge naturally from repeating periods of high and low refractive index films.On the other hand, our Fourier filters are based on a sinusoidally shaped transmission spectrum in Eq. ( 6) that rolls on and off smoothly.Although in principle any spectral shape can be achieved by incorporating a large number of layers and a numerical optimization method, the periodic nature of the spectrum suggests the use of a simple Fabry-Perot resonator as a basic element.This is the design approach taken in this work.Though periodic, Fabry-Perot resonators typically contain sharp peaks of transmission surrounded by broad high reflection bands, expressed by: where T a and T b are the transmission, and R a and R b are the termination reflectivities of each side of cavity (typically one side is air, the other side is the substrate), d s and n s are the physical thickness and refractive index of the spacer layer, and θ s is the oblique incident angle [11].
Assuming normal incidence, the termination reflectivities dictate the minimum transmission between the peak values, viz.they dictate the filter contrast χ which can be expressed by: A perfect sinusoidal transmission filter with same contrast χ can be written as where max T f p − χ 2 plays the role of α k in Eq. ( 6).See blue lines in Fig. 4. Comparing the Fabry-Perot filter to the perfect sinusoidal filter, we can define an error function ε to indicate the deviation of the Fabry-Perot filter from the desired sinusoidal function.In the 450 nm to 900 nm wavelength range, this error is represented by: Assuming silica glass as the substrate and air as the medium above the film, we can sweep the refractive index of the Fabry-Perot cavity to compute T a , T b , R a and R b and then plot ε as a function of the filter contrast χ.This is shown in Fig. 2. The transmission spectrum approaches a perfect sinusoid as the contrast drops to zero (which occurs when the refractive index of the cavity material approaches the substrate index).Therefore, there exists a trade-off between contrast and the sinusoidal shape.In addition, the optical dispersion of the cavity material also plays an important role in the periodicity of the transmission peaks.Ideally, a flat dispersion is required to get a perfect fit to a sinusoid across the entire spectral range of interest.In this work we chose ZnS as the cavity material.It is a fairly stable dielectric that is easy to deposit and has a relatively high refractive index of around 2.35 in the visible spectrum.The corresponding reflectivities are R 1 = 0.162 and R 2 = 0.048.The contrast χ is 28.8%, and the error can be found from Fig. 2 to be 1.6%.A different material such as TiO 2 can provide a larger contrast of 34% (using the refractive index at 650 nm), but it has a slightly greater dispersion of 12.34% between 450nm and 900nm, leading to a larger deviation from the sinusoidal function.Al 2 O 3 , on the other hand, has a much smaller dispersion of 1%, but its contrast will be significantly smaller at 8.4%.Therefore, ZnS is provides a reasonable compromise between performance and deviation from the ideal spectral shape.Finally, we are able to change the periodicity of the filters from 0.5 sinusoidal period to 3.0 sinusoidal periods across the spectral range of 450 nm to 900 nm by changing the thickness of cavity.Table 1 shows the design thickness of the ZnS layers.

Fabrication of Fourier spectral filters
All the filters were fabricated by depositing the ZnS dielectric thin films on fused silica substrates using RF magnetron sputtering.The uncoated substrates were cleaned and dried with nitrogen, and then placed in the sputtering chamber, which was then evacuated to a base pressure below 10 −6 Torr.A 3-minute pre-sputter step was used with the shutter closed to clean the ZnS target surface, and the deposition was performed with an RF discharge power of 100W in an Argon plasma with a working pressure of 8 mTorr.The film thicknesses were monitored on a quartz crystal thickness monitor and verified with a stylus profilometer after the deposition.
The six fabricated filters are shown in Fig. 3(b) mounted on a rotating filter wheel.The optical performance of the filters were measured using the Filmetrics F10-VC spectrometer.The transmission spectra of all six Fourier spectral filters are shown in Fig. 4 in wavenumber domain.The blue line is the pure sinusoidal transmission; green line is the design transmission; the red line is the measured transmission.We can verify that the measured spectral shapes match very well with the designs and adequately with the pure sinusoid, especially with respect to the positions of peaks and valleys.The small spectral shifts are attributable primarily to the material dispersion.

Experiments
A custom multispectral imaging test bed designed for operating ranges of 450nm to 900nm was built with our fabricated filters, as shown in Fig. 3.The system is comprised of three main components-a grayscale computer vision camera (OptiTrack V120 Slim); two cutoff filters (a high-pass filter at 450 nm and a low-pass filter at 900 nm); 6 sinusoidal filters mounted on rotating filter wheels.By time-multiplexing, we captured seven measurements (K = 7 in Eq. ( 6)) for each scene-one without filter (Z 0 (n n n)) and six filters ranging from 0.5 to 3.0 sinusoid periods.We also captured a blank image (image with no input) to determine the zero offset of the sensor.We also prepared a "target" comprised of narrowband LEDs peaking at 8 different wavelengths (460 nm, 515 nm, 565 nm, 590 nm, 627 nm, 655 nm, 850 nm and 880 nm), as shown in Fig. 3(c).Using this target, the experiments below verify the claim that spectral narrowband peaks can be recovered from our multispectral system based on sinusoidal filters.
In our first experiment, we measured the system response of LEDs individually.The spectra reconstructed from actual sensor measurements using Eq. ( 5) are shown in Fig. 5(a), and Table 2 shows the differences between estimated and the actual signal peaks.The peak error is small within 500-800nm range.Our results are consistent with the same experiment repeated with simulated sensor measurements (see Fig. 5

(b) and Table 2 last column).
The errors in the peak locations outside of the 500-800nm range can be explained by the results of our second experiment, where we captured a mixture of two LED lights using a  Note that 850nm emission is outside of visible range, so it does not appear in color rendering.
diffuser.When we combined 515nm and 850nm LEDs, the recovered spectrum in Fig. 7(a) clearly formed two peaks at the desired wavelengths.However, the combination of 565nm and 655nm LEDs yielded less distinct peaks-see red plot in Fig. 7(b).This is due to the fact that peaks at 565nm and 655nm broadened by the sinc function I 7 in (5) contaminate each other.
To illustrate this point further, suppose we make fewer measurements (K = 6) by ignoring Z 6 (n n n).Since sinc function I 6 is even broader, peaks at 565nm and 655nm merge and become indistinguishable.See the blue plot in Fig. 7(b).We conclude that the number of filters effectively determines the "spectral resolution" of the multispectral image sensor-I K becomes a sharper sinc function when K increases, making it possible to resolve multiple spectral features that are close in wavenumber or wavelength.This phenomenon explains the reasons for why the single peak experiments had larger errors in spectral ranges of 450-500nm and 800-900nm in Table 2. Recall that the reduced measurement strategy of Eq. ( 3) is only valid between σ ∈ (900nm −1 , 450nm −1 ).By Fourier periodicity, X(n n n, σ ) in Eq. ( 3) is a symmetrically extended-in the sense that σ = 460nm −1 , 850nm −1 , 880nm −1 peaks are repeated at σ = 440nm −1 , 956nm −1 , 921nm −1 , respectively (440nm −1 = 2 × 450nm −1 − 460nm −1 , etc.).Since 440nm and 460nm peaks are very close, the broadened peaks contaminate each other (same with 850nm and 956nm; or 880nm and 921nm).To overcome this problem, one can increase the number of filters.Alternatively, one can decrease the OPD sample rate ∆ξ , though this comes at the sacrifice of the sensor operating range.One drawback to our approach is Gibbs phenomenon, or extra ripples stemming from the oscillations of the sinusoids.This is only a problem when the spectral features are sharp (such as the narrowband specta of LED's)-for example, there are false peaks near 650nm in Fig. 7(a) and at 450nm in Fig. 7(b) and 7(c).When spectrum is represented in wavelength ( X(n n n, λ −1 )λ −2 ) instead of in wavenumber (σ = 1/λ ), a small Gibbs ripple near 450nm is likely to appear bigger, while 900nm appear smaller (as evidenced in Fig. 5).Gibbs phenomenon is reduced when the number of filters increase or the OPD sample rate ∆ξ decrease.Finally, the above analysis above generalizes to more complex and/or broadband spectra.Since spectrum reconstruction via Fourier transform in Eq. ( 5) is completely linear, reconstruction of any spectral scene is a linear combination of system responses at single wavenumber responses X(n n n, σ ) integrated over σ ∈ (900nm −1 , 450nm −1 ).For example, reconstruction from 3 LEDs combined in Fig. 7(c) are simply a linear combination of the single LED spectra in Fig. 5(a) combined together.

Conclusions
In this paper, we have developed and demonstrated the concept of Fourier multispectral imaging.We designed and built periodic filters to approximate the sinusoidal transmittance.Using these filters, we made finite OPD measurements, and the spectra were reconstructed by the forward Fourier transform.We theoretically proved and empirically verified that finite OPD samples preserve the spectral features such as transmission and absorption peaks.This makes the Fourier multispectral imaging technique an attractive alternative to narrowband filter-based multispectral imaging.

Fig. 2 .
Fig. 2. Percent error between the Fabry-Perot filter function and the ideal sinusoidal function as a function of the filter contrast χ.

Fig. 4 .Fig. 5 .
Fig. 4. The blue lines show the desired pure sinusoidal transmittance function.The green lines are the designed filter spectra that best approximates the blue line.The red lines are measured transmission spectra from the fabricated filters.

Fig. 6 .
Fig. 6.Raw sensor measurements under different filters and reconstructed image when multiple LEDs are lit.LED peaks (top row) at 515nm and 850nm, (middle row) at 565nm and 655nm, and (bottom row) at 565nm, 655nm, and 850nm.(a-g) The seven multispectral measurements {Z 0 , . . ., Z 6 } taken.(h) Color image rendered from reconstructed spectra.Note that 850nm emission is outside of visible range, so it does not appear in color rendering.

Fig. 7 .
Fig.7.Reconstructed spectra in Fig.6with the specified LEDs turned on.The blue and red lines represent reconstructions from 6 and 7 measurements, respectively.In (b-c) 565nm and 655nm LED can be distinguished using K = 7 measurements, but peaks are indistinguishable with only K = 6 measurements.The ability to resolve multiple spectral features that are close in wavenumber or wavelength determines the number of filters needed.

Table 1 .
Designed ZnS thickness of Fourier spectral filters

Table 2 .
Comparisons of reconstructed and actual spectral peaks of LEDs in Fig.5.