Image Reconstruction from SPDs
Images can be reconstructed from SPDs standard methods. [7] The preceding discussion has neglected the subject of units/scaling due to the associative nature of the equations. However, scaling becomes important in image reconstruction due to the fixed limits of the color space.
Any set of observer functions (e.g. for a population, camera, individual, etc.) can be used for image reconstruction. The most common of these is the 1931 CIE 2° XYZ observer. [8] Taking the dot product of the SPD at each pixel with each of the observer sensitivity curves, an XYZ tristimulus color is obtained. These XYZ colors are then normalized 0-1 and converted to RGB for display.
Standard tristimulus observer functions for image reconstruction from SPDs
Because SPDs are relative, not absolute, the sensor has no ability to set the black point and white point in the image. Instead, the operator may manually specify where these limits occur according to their knowledge of the scene, or aesthetic goals for the image. They may also elect to scale their image to the limits of the color space for maximum contrast, although this is generally unlikely to produce perceptually accurate colors.
A simple default scaling method is suggested here:
- Convert the image to grayscale.
- Calculate the probability density function (PDF) of the image pixel values.
- Calculate the cumulative distribution function (CDF) as the cumulative sum of the PDF.
- Look up the values corresponding to 1% and 99% in the CDF
- Using imadjust or equivalent, scale the value range of the image so that 1% and 99% on the CDF correspond to 5% and 95% on the value range
These values are good for default general-purpose scaling, but may be adjusted depending on the specific scene being imaged.
Computational Approach
All computation described in this article is implemented in MATLAB. Although SPDs may be calculated pixel-by-pixel, this approach scales poorly to multiple-megapixel images. Instead, SPDs are calculated at the datacube level using reshape and repmat to eliminate multiple for loops and leverage MATLAB's matrix math efficiency. This works as follows: - Initialize an SPD datacube with zeros for all entries, measuring m-by-n-w, for photos of resolution n-by-m and w wavelengths
- Pick an ST curve
- Retrieve the color channel of the photo corresponding to this ST curve
- Reshape the ST curve along the 3rd (wavelength) dimension with reshape
- Repeat the ST curve at each pixel with repmat
- Repeat the color channel at each wavelength with repmat
- Calculate the element-wise product of the color datacube and the ST datacube
- Add this result to the SPD datacube
- Repeat back to Step 2 for all remaining ST curves
- Divide the SPD datacube by the sum of all ST curves, similarly reshaped and repeated
- Multiply the SPD datacube by the illuminant and calibration curves, similarly reshaped and repeated
Depending on the camera and filter set, some wavelengths may have zero or near-zero sensitivity. This commonly occurs at the fringes of the visible light spectrum, where the camera is least sensitive. At these wavelengths, the SPD is undefined, or nearly so. To address this, wavelength channels corresponding to a sensitivity of < 1% of the peak value are discarded, and set equal to the closest non-discarded wavelength channel. Depending on the specific camera-filter system, it may be preferable in other cases to interpolate or extrapolate instead.
This "insensitivity" check is performed on U, the sum of the ST curves over the st index, shown below. It may be thought of as the overall sensitivity of the camera/filter system, as a function of wavelength. For a truly non-overlapping filter set, e.g. for a narrow bandpass filter set, U would have a sawtooth shape, and interpolation would be preferable in the insensitive regions between peaks. For the regular bandpass filter set under consideration, overlap occurs, and insensitivity occurs only at the extremes of the domain, where the SPDs from the closest "sensitive" wavelength are simply copied, or held.
Overall sensitivity, used to normalize and determine of where to interpolate/extrapolate
Filter Set Selection
As discussed in Seeoung et al. [4], real-world reflectance spectra can empirically be decomposed into approximately 8 principal components, supposing > 99% of variance is to be explained. SPDs differ from reflectance by a factor of the illuminant, which by inspection is of lower dimensionality for everyday cases. [8] [9]
From the perspective of sensor design, this indicates that SPDs should be well-approximated by 8 values distributed evenly over the visible wavelength band of roughly 400 to 700 nm, equivalent to a spectral resolution of approximately 43 nm. Between these wavelength values, interpolation can be used to good approximation. For an image sensor with monochromatic sensitivity, this would necessitate at least 8 filters. However, the trichromate sensitivity of ordinary cameras can be used to decrease the filter count.
The principal components of a camera/filter set system can be evaluated using the ST matrix shown above:
which expands as follows:
where 1, 2, 3, ... λ denotes the wavelength index and a, b, c, ... st denotes the ST index.
At an arbitrary pixel, the ST matrix can be used to set up an over-defined system of linear equations per the photo capture expression:
where COL is a vector of color measurements, and SPD is a vector of the spectral power distribution.
This system of linear equations is equivalent to the canonical Ax=b formulation with a known measurement vector b (COL), known transformation matrix A (ST), and unknown predictor vector x (SPD). It is mathematically possible to solve this over-defined system directly, e.g. with a linear least-squares (LLS) method:
However, the ST matrix tends to be poorly conditioned, producing non-physical and non-convergent results for SPD if using LLS. LLS is also prohibitively slow at megapixel scale.
Nevertheless, the ST matrix is still useful for principal component analysis (PCA), to estimate the spectral resolution of the camera / filter set system. Because PCA is a function only of the transformation matrix ST, it can be performed in the absence of any photo data, and is a powerful tool for evaluating filter sets.
It can be shown using PCA that the ideal filter set has transmission spectra with minimal overlap and maximum uniformity. If the spectra do not overlap, the filters are orthogonal in SPD vector space, and their utility in terms of PC per filter is maximized. It follows that ideal filters have "spike" spectra that measure a single wavelength only. In practice, this is best realized with a series of bandpass filters with minimal FWHM.
Using this PCA approach, several filter sets were evaluated, ranging in quantity from 5 to 12, and in cost from $30 to $3,500. These included options from Thorlabs, Edmund Optics, MidOpt, and PixelTeq. These fell roughly into the following tiers:
Comparison of filter set tiers by key cost and performance metrics; selection highlighted
The tiers differ significantly in terms of cost and performance. Gel-like color filters are very cheap, [2] but because their transmission spectra are arbitrary and usually non-zero nowhere, they overlap significantly, and it is practically difficult to develop more than 5 PC. Further, they usually lack published transmission data. Bandpass filters with 60-90 FWHM are more expensive but get much closer to the target of 8 PC. Their relatively broad FWHM limits the maximum number of useful filters than can be used to ~5, and the PC count to ~7. Narrow bandpass filters with 20-40 FWHM are more expensive yet but offer greater optionality. They can be made to develop as many as 14 PC. 8 PC is also achievable, but still at substantially greater cost than with regular bandpass filters and 7 PC, at ~$1,400 vs.$500.
Size and form factor/deployability were also considered. The ~43 mm diameter size required is relatively large relative to typical offerings, which commonly are limited to 1" [25.4 mm]. Many filter kits come in large, bulky cases, or are unmounted.
Ultimately, the MidOpt FS100 was selected. At $500 and 10 filters, it is an outlier in value, despite only 5 of the 10 filters being useful for visible light imaging. The swiveling fan form factor is ideal for no-contact lens filtering, as it can be operated with one hand, protects the filters, and is extremely compact when collapsed.
Several filters in the FS100 kit are not useful for this application as they either operate outside the visible light range, are polarizing filters, or excessively overlap another filter. The useful subset of filters were removed and re-bound using a shorter brass binding post:
Subset of useful filters from MidOpt FS100
The filter colors appear artificially warm in this picture due to their red-orange reflective coatings. The final filter subset is:
- MidOpt BP470, FWHM 85 nm, blue
- MidOpt BP525, FWHM 80 nm, light green
- MidOpt BP590, FWHM 75 nm, orange
- MidOpt BP635, FWHM 60 nm, light red
- MidOpt BP660, FWHM 65 nm, dark red
Bandpass transmission spectra for filters chosen from MidOpt FS100 kit
Colors for filters chosen from MidOpt FS100 kit - for visual check only, not used in calculation
This filter set is shown to have 7 PC of diminishing explanative power, using the same standard as Seoung et al.: [4]
Principal component analysis of camera/filter system, showing ~7 PC
Photographic Aspects
Camera Selection
It is a prerequisite for this method that the spectral sensitivity of the camera is known. This data is not typically reported by the manufacturer, nor can it be measured without specialized equipment. Fortunately, literature provides data for some common models, most prominently Jiang et al. who published a database of 28 models in 2013, which is used in this study. [1]
In some cases, it may be acceptable to use a similar camera's data, particularly if the manufacturer is the same and the model year/type is similar. Because camera spectral sensitivities are designed to mimic the human eye, they generally resemble standard observer functions with minor deviations.
In this study, photos are captured with a Canon 650D, with sensitivity data given by the similar Canon 600D.
Photo Capture
Minimizing perturbation of the camera while capturing the photo stack is critical. To this end, a tripod and remote shutter are used. The filters are neither attached to nor in contact with the camera at any time; they are simply held in front of the lens by hand at a distance of approximately 1/4". A prime (i.e. non-zoom, or fixed focal length) lens is preferable as it has one fewer degree of freedom than a zoom lens, and is less prone to perturbation and backlash. Prime lenses also have smaller apertures than zoom lenses, which enables use of smaller filters. A Canon 40 mm prime lens is used in this study, chosen for its generality and naturalistic field of view. Accounting for the crop-sensor of the Canon 650D, the effective focal length is approximately 65 mm.
Photography setup: filter fan in left hand, remote shutter in right hand; no contact with camera
Since this method deliberately avoids the RAW format, care must be taken to avoid saturating the limits of the color space while capturing photos, particularly by overexposure. Saturation occurs in color channels individually before the corresponding grayscale image saturates; therefore the R, G, and B histograms should be inspected individually, rather than the grayscale histogram by itself.
ISO and exposure time should be minimized to reduce sensor noise and blurriness, while still ensuring that the image adequately spans the color space. The specific settings for ISO and exposure time will vary depending on the subject and imaging goals. White balance should be set to match the scene as best as possible.
It is advisable to take several test photos without a filter to correct the camera settings before capturing the photo stack. The goal should be to get the unfiltered photo just shy of overexposure, so that the filtered photos have maximum dynamic range without saturation. Because the MidOpt filters have a very high peak transmission of 90-100%, no adjustment of ISO or shutter speed is necessary between the test photos and the stack photos. These test photos also serve as a qualitative check for the reconstructed images.
Validation and Calibration
SPD estimates are compared against samples of known properties - specifically, a color chart consisting of 108 samples of oil paint mixtures spanning the color gamut, originally developed by painter Richard Schmid. [10] [11] SPDs are estimated by capturing a filtered photo stack of the Schmid chart under noon daylight, corresponding to the CIE D65 illuminant. SPD ground truth is calculated as the product the D65 illuminant and the spectral reflectance of each swatch, as measured by the Spectro 1 spectrophotometer. The quality of SPD estimation is assessed by comparing the estimated and measured SPDs for each swatch.
Comparison of measured vs. estimated SPD at each swatch for un-calibrated system
From this comparison, a calibration curve is developed. At each swatch, the ratio of estimated to measured SPD is calculated, i.e. the required gain to eliminate error between measurement and estimation. These gains are averaged at each wavelength for all swatches, resulting in the following calibration curve:
Overall error signal for Schmid chart before calibration
To maintain generality and avoid overfitting, a simple linear calibration curve is fit to this curve:
Fitted calibration curve to minimize error signal for Schmid chart
The Schmid chart is then reevaluated using the calibration curve:
Comparison of measured vs. estimated SPD at each swatch for calibrated system
The fit between measurement and estimation is much improved. The average error is centered on 1.00, indicating that the linear fit has the intended effect.
Overall error signal for Schmid chart after calibration
Results
Several results are shown. All are calculated with a D65 illuminant, the Schmid calibration curve applied, a wavelength resolution of 20 nm, and value scaled to 5% and 95% with a nominal gamma (linearity) value of 1.00. For each, the unfiltered/test (left) and reconstructed (right) images are compared, and a mesh of characteristic spectra are shown.
Schmid Color Chart
Cabbage
John W. Weeks Bridge
Little Fresh Pond 1
Sugar Maple
Little Fresh Pond 2
John W. Weeks Bridge Entrance
Belmont Victory Gardens, Rock Meadow
JFK Memorial Park 1
JFK Memorial Park 2
Desk
Conclusions
A simple method for hyperspectral imaging with a commodity camera is demonstrated using only commercially-available hardware. Accounting for the camera's Bayer filter reduces the number of lens filters required by increasing the number of principal components in the photo capture transformation matrix, ST. For the specific choice of a Canon 650D camera and the MidOpt FS100 filter kit, reasonable agreement is shown for samples of known properties between measured vs. estimated SPDs. Future work may examine in greater depth the filter selection sub-problem to improve SPD estimation and/or reduce cost, and/or study differences between unfiltered and reconstructed images.
Source Code
A MATLAB GUI was developed to facilitate rapid adjustment and iteration of SPD estimation and image reconstruction parameters. It currently supports 2 filter sets, 7 illuminants, 28 cameras, and 2 observers, and can readily be expanded as needed. Full source code may be found on GitHub.
MATLAB GUI for creation of SPDs and reconstructed images
References
- Jiang et al. What is the Space of Spectral Sensitivity Functions for Digital Color Cameras? IEEE Workshop on the Applications of Computer Vision (WACV), 2013.
- Leeuw et al. In situ Measurements of Phytoplankton Fluorescence Using Low Cost Electronics. Sensors 2013, 13, 7872-7883.
- Antonino Cosentino. Multispectral Imaging of Pigments with a Digital Camera and 12 Interferential Filters. e-Preservation Science, 2015.
- Seoung et al. Do It Yourself Hyperspectral Imaging with Everyday Digital Cameras. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016.
- Habel et al. Practical Spectral Photography. Eurographics, Volume 31, Number 2, 2012.
- Baek et al. Compact Single-Shot Hyperspectral Imaging Using a Prism. ACM Transactions on Graphics, Vol. 36, No. 6, Article 217. Publication date: November 2017.
- Lindbloom. Spectral Computation of XYZ. Revised Sat, 08 Apr 2017 18:57:27 GMT.
- Colour & Vision Research Laboratory, University College London. Standard Illuminants.
- PublicResource.org. CIE 15: Supplementary Tables [10 CFR 430 Subpart B, App. R, 4.1.1], cie.15.2004.tables.xls.
- Dichter. How Paints Mix. 2019.
- Schmid. Alla Prima II. Page 219. Fourth Printing, March 2018.