Optimized spectral filter design enables more accurate estimation of oxygen saturation in spectral imaging

Oxygen saturation (SO2) in tissue is a crucially important physiological parameter with ubiquitous clinical utility in diagnosis, treatment, and monitoring, as well as widespread use as an invaluable preclinical research tool. Multispectral imaging can be used to visualize SO2 non-invasively, non-destructively and without contact in real-time using narrow spectral filter sets, but typically, these spectral filter sets are poorly suited to a specific clinical task, application, or tissue type. In this work, we demonstrate the merit of optimizing spectral filter sets for more accurate estimation of SO2. Using tissue modelling and simulated multispectral imaging, we demonstrate filter optimization reduces the root-mean-square-error (RMSE) in estimating SO2 by up to 37% compared with evenly spaced filters. Moreover, we demonstrate up to a 79% decrease in RMSE for optimized filter sets compared with filter sets chosen to minimize mutual information. Wider adoption of this approach will result in more effective multispectral imaging systems that can address specific clinical needs and consequently, more widespread adoption of multispectral imaging technologies in disease diagnosis and treatment.


Introduction
Oxygen saturation (SO 2 ) is a crucially important physiological parameter with ubiquitous clinical utility in diagnosis, treatment, and monitoring, as well as widespread use as an invaluable preclinical research tool. Typically, a significant decrease in SO 2 indicates a disruption of normal biological function. Traditionally, SO 2 is measured non-invasively using pulse oximetry, or invasively using bedside equipment, such as spectrophotometers, that measure an extracted blood sample. The latter requires a relatively large blood sample, limiting its applicability to measuring arterial or venous blood, and carries inherent drawbacks associated with invasive procedures including risk of infection and use of expensive single-use equipment in addition to workflow challenges arising from the sample acquisition. Neither approach provides spatially resolved information and they are thus unable to resolve local variations in blood saturation that might be useful in assessing certain pathologies. In recent years, research into real-time non-invasive, non-contact optical techniques for visualizing SO 2 has gained momentum and these techniques have found application in a wide range of indications.
Intraoperative visualization of SO 2 is valuable across a wide range of surgical specialties. Insufficient perfusion of tissues can result in ischemic injury, reduced viability of repaired tissue and poor healing. Some diseased tissue has a different vascular profile to its healthy form and further, the identification of blood vessels, which may be in unexpected locations due to disease or anatomical variations, is critical in avoiding accidental injuries that might otherwise result in prolonged procedures and serious or life-threatening complications [1]. Moreover, the success of anastomosis, the surgical attachment of two luminal structures [2][3][4], and the creation of skin flaps [5,6], one of the most common surgical techniques for repairing missing or damaged tissues, rely on ensuring proper tissue perfusion. In the brain, dynamic monitoring of SO 2 is particularly important to safe and effective cerebrovascular reconstruction [7].
Monitoring SO 2 is also useful for transplant surgery: for monitoring donor organ perfusion during warm ischemia, whether this be in donor, recipient or during normothermic machine perfusion [8]; and for monitoring reperfusion during surgery [9], ensuring surgical attachment of the organ and its associated vasculature is successful, and promoting long-term viability of the graft.
Because oxygen supply is a fundamental factor for healing, non-invasive visualization of SO 2 also plays a key role in objective assessment of wound damage and healing potential [10] for heat burns [11,12], chemical burns [13] and radiation burns [14], as well as for diagnosis and monitoring of chronic wounds [15], such as chronic skin ulcers [16] and diabetic foot ulcers [17,18], and for monitoring hemodynamic disorders such as scleroderma and Dupuytren's contracture [19].
Imaging of SO 2 is particularly suited to cancer imaging as angiogenesis, the development of tumor-associated neovasculature, is one of the hallmarks of disease [20]. Mapping of SO 2 might exploit these changes for early detection of cancer, defining lesions based on changes in tissue oxygenation or vascular properties [21]. It can also be used to characterize hypoxia in solid tumors, which is related to treatment-resistance and selection of appropriate therapeutic strategies [22]. Similarly, it has also been used to aid diagnosis and monitoring of pathologic conditions of the retina and optic nerve, where loss of normal oxygen supply is believed to play an important role in disease [23,24].
Beyond application in the clinical setting, visualization of SO 2 is useful in basic research, for example, mapping of hemodynamic response in the brain to understand brain organization and processing [25] or response to blast-induced traumatic brain injury [26], and for assessment and monitoring in tissue engineering.
Hyperspectral imaging [27][28][29] (HSI) has the potential to non-invasively measure SO 2 based on the distinct absorption and scattering spectra of oxyhemoglobin (HbO 2 ) and deoxyhemoglobin (Hb). HSI captures both spatial (x,y) and spectral (wavelength,λ) information to acquire an image (hyper)cube (x,y,λ). Spectral unmixing is then used to estimate the abundances of HbO 2 and Hb in each image pixel based on the captured spectra. By measuring their relative abundance, SO 2 can be estimated, while the sum of abundances indicates the total blood volume. This approach benefits from being real-time, non-invasive, non-contact and label-free.
Four main acquisition approaches are used: point scanning, where complete wavelength information is captured from each image point sequentially, line-scanning, where complete wavelength information is acquired from each image line sequentially, wavelength scanning, where an entire image is captured at each wavelength sequentially, and snapshot approaches, where the entire image cube is captured in a single snapshot. With all approaches, there is a trade-off between speed, spatial resolution and spectral resolution. Devices with high spectral and spatial resolution tend to be slow, bulky, costly and susceptible to misalignment -unsuitable for a clinical environment where low-cost, robust, real-time imaging is desirable [30]. Additionally, capturing more wavelengths at higher resolution leads to larger image cubes, requiring larger storage systems, longer save times, and slower classification and display. For these reasons, techniques intended for clinical translation tend to capture a reduced number of narrow wavelength 'bands' and thus perform 'multispectral imaging' (though there is no universally agreed limit on the number of bands required to distinguish hyper/multi-spectral).
Multispectral imaging is typically implemented by wavelength scanning, either using a fast filter wheel or tunable light source, or more recently by using spectrally resolved detector arrays (SRDAs), exploiting spectral filters deposited directly onto the imaging detector in a mosaic pattern. Whichever acquisition regime is used, a fundamental question remains: which wavelengths should be captured to visualize SO 2 with maximum accuracy?
Typically, this question has been avoided in the development of biomedical multispectral imaging, with developers deploying general purpose 'off-the-shelf' sensors with ∼10 evenly spaced narrow bands across the wavelength range of interest [27][28][29]. One of the key challenges is the lack of reliable gold-standard datasets of biological tissue spectra [30], so those developers that have attempted to address the question of spectral band selection have done so using a 2-stage development process, first capturing high resolution data from their tissue of interest using slow scanning HSI devices, before using the high-resolution data to perform spectral band selection for MSI.
This approach has seen some promising successes. Wirkert et al. (2014) analyzed surgical image cubes comprised of 30 wavelengths and used an information-theory-based approach to identify seven optimal bands for SO 2 estimation, allowing them to perform MSI with a fast filter-wheel device [31]. Waterhouse et al. (2021) analyzed diffuse reflectance spectra to identify three bands that would optimally display contrast between Barrett's esophagus and cancer in the esophagus, achieving >12-fold contrast enhancement [21]. Perhaps the biggest success of the 2-stage approach, and certainly the most clinically used in routine care, is demonstrated by the invention of narrow band imaging (NBI). The 2 filters used for NBI were selected from 9 off-the-shelf filters for their ability to enhance contrast for vasculature in the human tongue [32]. Even with this relatively crude selection of only two bands, NBI has proven advantageous in the detection and characterization of early Barrett's esophagus-related neoplasia. It was the first advanced imaging technique to meet the requirements for recommendation in Barrett's esophagus surveillance [33] and has been successfully translated into routine clinical practice.
From a statistical point of view, spectral band selection is viewed as an optimization problem in which we seek a subset of spectral bands that capture most of the information for a particular hyperspectral signal, whilst removing noise and spectral redundancy (highly correlated wavelengths). However, for clinical application, it may not be necessary to collect all the information, as much of it may be non-discriminatory. In this sense, spectral band optimization for clinical application aims to eliminate noise and spectral redundancy whilst preserving discriminant (or diagnostic) information only.
Spectral band selection has been developed for decades in the signal processing community [34], particularly for application in land surveillance, but few of the techniques that have been developed in this field have been applied to biomedical imaging problems. Where spectral band optimization has been applied to biomedical problems, researchers aimed to maximize accurate classification of tissues. However, rather than directly maximize classification performance itself, which is computationally expensive as it requires classifiers be trained and tested for each filter set, the studies used indirect 'functions of merit'. Some studies maximized mutual information [31,35,36], which neglects the above-mentioned subtlety that not all information is useful information; others maximized differences between the characteristic spectra of classes, such as root-mean-square-difference or variance [37], but this is not applicable to continuous gradient problems such as estimating SO 2 . Borrowing vocabulary from feature selection for machine learning, both methods are 'filter' methods rather than 'wrapper' methods, meaning that whilst each selected band contributes to maximizing the chosen function of merit, there is no guarantee that selected filters will result in the best performance once applied to the final classification problem.
In this paper, we compare 'off-the-shelf' filter sets to those selected using optimization approaches; both indirect optimization by maximization of root-mean-square-difference, spectral angle and mutual information, and direct optimization by minimization of error in estimated abundance. Through these comparisons, we demonstrate the merits of direct application-specific spectral band selection and advocate for future work in this area to push multispectral imaging to reach its full potential to impact clinical care through imaging systems tailored to the clinical target and needs.

Framework for optimizing spectral filter sets
An overview of our approach to optimizing spectral filter sets is shown in Fig. 1. Briefly, a tissue signal hypercube, S(x,y,λ), is modelled using a ground truth abundance map, A true , to determine the abundance of oxy-and deoxy-hemoglobin in each pixel, and an empirical model to calculate the spectrum of diffusely reflected light in each pixel using the absorption and scattering coefficients of oxy-and deoxy-hemoglobin, µ a and µ' s , and the spectrum of the light source, L(λ). Multispectral images, I(x,y), are simulated by propagating the tissue hypercube through n filters defined by center wavelengths λ c = [λ c,1 , . . . , λ c,n ] and full width half maxima w = [w 1 , . . . , w n ]. Pixel-wise spectral unmixing of the simulated images results in an estimated abundance map, A est. (x,y) which is compared to the ground truth abundance map to determine the root-mean-square-error (RMSE). The RMSE is minimized by adjusting filter parameters in an optimization loop until predefined stopping criteria are met, resulting in an optimized filter set.

Fig. 1. Overview of the spectral filter optimization framework.
A tissue signal hypercube, S(x,y,λ), is modelled using a ground truth abundance map, A true , to determine the abundance of oxy-and deoxy-hemoglobin in each pixel, and an empirical model to calculate the spectrum of diffusely reflected light in each pixel using the absorption and scattering coefficients of oxy-and deoxy-hemoglobin, µ a and µ' s , and the spectrum of the light source, L(λ). Multispectral images, I(x,y), are simulated by propagating the tissue hypercube through n filters defined by center wavelengths λ c = [λ c,1 , . . . , λ c,n ] and full width half maxima w = [w 1 , . . . , w n ]. Pixel-wise spectral unmixing of the simulated images results in an estimated abundance map, A est. (x,y) and this is compared to the ground truth abundance map to determine the root-mean-square-error (RMSE). The RMSE is minimized by adjusting the filter parameters in an optimization loop until predefined stopping criteria are met, resulting in an optimized filter set.

Tissue modelling to generate signal hypercubes
Signal hypercubes were simulated as follows. A ground truth abundance map, A true (x,y,i), is defined such that, where x and y are spatial coordinates in the map, A true HbO 2 (x, y) are the ground truth abundances of oxyhemoglobin (HbO 2 ) and deoxyhemoglobin (Hb) at the point (x,y) respectively. This map had a size of 21 × 21 pixels. The oxygen saturation, SO true 2 (x, y), at the point (x,y) is determined as, and the total abundance of hemoglobin (THb) is defined as, The map was defined such that a linear gradient of SO 2 (from 0 to 1) was present horizontally across the map, while a linear gradient of THb (from 0 to 1) was present vertically across the map, thus ensuring all combinations of SO 2 and THb are evenly represented in the map.
To simulate the signal detected from diffuse reflectance imaging of the ground truth abundance map, we considered several methods of calculating tissue reflectance from absorption and scattering coefficients. These included calculation from the effective attenuation coefficient [38], a model of diffuse reflectance under diffuse illumination using three flux theory [38], a model of diffuse reflectance using a broad beam model [38] and an empirical model [39]. By comparing the modelled reflectance to real tissue spectra [21], we determined the empirical model to be the best approximation for diffuse reflectance imaging of in vivo biological tissue (Fig. S1).
For each pixel in the abundance map, the empirical model is used to calculate the signal for each pixel of a signal hypercube, where L(λ) [ [40]. The values of k 1 = 0.26 and k 2 = 14 were determined by fitting the empirical model to the mean of 320 spectra captured from in vivo human esophageal tissue in a previous study [21] (Fig. S1). Gaussian white noise was added to each spectrum using the MATLAB function 'awgn' with a signal to noise ratio of 100. This results in a hypercube of spectral data, S(x,y,λ).

Simulation of multispectral imaging
For simulated imaging, a set of i = 1, . . . , n spectral filters are defined as, where, where N is a normalization factor that ensures the area under the curve is equal to 1, λ c,i is the center wavelength of the ith filter, and w i is the full width half maximum (FWHM) of the ith filter such that, are the center wavelengths and FWHMs of the filter set F. To simulate imaging with the filter set, the signal hypercube, S(x,y,λ), an image cube is calculated, where, Zero-mean gaussian white noise was added to each image using the MATLAB function 'imnoise' with variance 5 × 10 −5 . To simulate the auto exposure function present in most imaging systems, the image cube was normalized to the maximum pixel value in the whole image cube.

Spectral unmixing to estimate oxygen saturation
For spectral unmixing, the reduced scattering coefficients, the absorption coefficients and the illumination spectrum are first propagated through the spectral filters to generate endmembers, While the ground truth absorption and scattering spectra might not be known for some applications, it is reasonable to suggest these could be measured in a calibration step, or otherwise calculated using databook values. For each pixel of the image cube, I(x,y), the estimated abundances of HbO 2 and Hb, A est. HbO 2 (x, y) and A est.
Hb (x, y), are estimated by least squares fitting of the image cube spectra with the empirical model described in Eq. (5-7), using the endmember spectra described in Eq. (14)(15)(16)(17)(18) as inputs. In other words, by minimization of the sum of square errors cost function, the abundances A est. HbO 2 (x, y) and A est .  Hb (x, y), can be estimated. The fitted scalar c accounts for normalization of the image cube. The fitted abundances form the estimated abundance map,

Merit function for performance calculation
To quantify the performance of a particular filter set, the root-mean-square-difference between the estimated abundances and the ground truth abundances was calculated as, where X and Y are the length and width of the image in number of pixels.

Optimization
To find the optimum filter set the merit function [Eq. (22)] is minimized. We consider the selection of n filters from a discrete set of filters with center wavelengths and FWHMs, Thus, an exhaustive approach to optimizing the filter set is not feasible.

Gradient descent
A gradient descent algorithm was used to minimize the merit function and thus find the optimum filter set. Gradient descent algorithms start from an initial estimate for the optimum parameters, in this case the filter properties λ c and w, typically obtained using a linear solver. From here, the algorithm takes repeated steps down the steepest local gradient in the merit function, with the aim of reaching a global minimum. The algorithm started with n filters with evenly spaced center wavelengths and w = 5 nm. The gradient descent algorithm moves each filter, i, sequentially. Briefly, the merit function [Eq. (22)] is calculated in a small region around the current filter position defined by, The filter is moved to the position in this region where the merit function is minimized. The occurs sequentially for each filter, i, until convergence (defined as no change in merit function for all filters, or a loop returning to a previously tested filter set). Throughout this process, overlapping filters are not allowed, 2.6.2. Genetic algorithm Gradient descent algorithms are susceptible to getting stuck in local minima and thus not finding a global minimum (Fig. 2). To overcome this limitation, a genetic algorithm was used for optimization. Inspired by biological evolution by natural selection, genetic algorithms start with an initial population of individuals, in this case a population of 20n (up to max 100) filter sets with filter parameters λ c and w generated uniformly at random within the bounds. For these, the 'fitness' is calculated using the merit function described in Eq. (22). Subsequent populations are generated based on the current population through three processes: mutation, introducing random changes in properties; crossover, combining the properties of pairs of parents; and automatic unchanged survival of the 'elite' members of the population. The algorithm terminates when the average change in the fitness value is below a set tolerance, indicating convergence. Overlap of bands was prevented by passing the inequality in Eq. (28) to the optimization function, effectively preventing such filter sets from being allowed in the population. The crossover fraction was set to 20% for the first 25 generations to search primarily at random via mutation, thus searching the entire optimization space, and 80% thereafter to search primarily via combination of the remaining members of the population.

Alternative merit functions based on maximizing the difference between HbO 2 and Hb
Previous attempts to optimize spectral filter sets have focused on maximizing the difference between the endmember spectra, such as root-mean-square-difference. These methods presume that maximizing the difference between the extremes of oxygenation (100% and 0%) will lead to greater performance in unmixing the signal contributions of HbO 2 and Hb, thus leading to more accurate estimation of SO 2 , but it is not clear if this is the case. Thus, to compare our approach of directly optimizing RMSE in SO 2 [Eq. (22)] to indirect methods, optimization was performed using two alternative performance metrics. Images were simulated as described in Sections 2.1-2.4, but with the SO 2 -gradient ground truth abundance map replaced with a binary ground truth abundance map consisting of one half having SO 2 = 1 and the other half having SO 2 = 0. This allowed the calculation of HbO 2 and Hb endmember spectra, HbO 2 and Hb, the mean endmember spectra in the SO 2 = 1 and SO 2 = 0 regions respectively. Optimization was performed as described in Section 2.6.2 using a genetic algorithm, with the performance metric in Section 2.5 [Eq. (22)] replaced by two alternative metrics: a metric based on the root-mean-square-difference (RMSD) between HbO 2 and Hb endmember spectra, and a metric based on the spectral angle (SA) between HbO 2 and Hb endmember spectra, where || || represents the Euclidean norm and | | represents the 1-norm.

Optimizing filter sets using mutual information
An alternative approach to selecting filters is to minimize the mutual information measured by the chosen filters. Briefly, for the signal hypercube, S, generated according to Eq. (5) with λ = [470, 485, . . . , 850] nm, a distance measure based on normalized mutual information [41] was defined as, with the normalized mutual information defined as,

NMI(S(λ i ), S(λ j )) = 2 H(S(λ i )) + H(S(λ j )) − H(S(λ i ), S(λ j )) H(S(λ i )) + H(S(λ j ))
where S(λ i ) and S(λ j ) are images in the signal hypercube at two different wavelengths λ i and λ j , H(S) are the marginal entropies of the images; and H(S i , S j ) is the joint entropy of the images S(λ i ) and S(λ j ). Entropies were calculated using the histogram method [41] with 100 bins. Based on this distance measure, agglomerative hierarchical clustering of the filter images was performed to generate n clusters, C. The resulting clusters represent mutually exclusive clusters of highly correlated filter images. To calculate distances between pairs of clusters, the average distance between all pairs of images in the two clusters was used.
Finally, the weight of each filter image within a cluster C is defined as, where n c is the number of filter images within the cluster and ε = 10 −6 to avoid singular values.
The filter corresponding to the filter image with the highest weight in the cluster is defined as the selected filter from the cluster. Thus, from n clusters, n filters are selected.

Assessment of optimized filter sets using test hypercubes
To compare filter sets, a series of simulated test hypercubes were prepared using ground truth maps according to Eq. (5-7) (Fig. 3). For biomedical applications, accurate visualization of oxygen saturation is crucial to enabling the detection of ischemia. To quantify the accuracy of each filter set in determining the oxygenation at different total hemoglobin abundances, a gradient map with a linear gradient of oxygen saturation, SO true 2 , (from 0 to 1) horizontally across the map, and a linear gradient of total hemoglobin abundance, A true THb , (from 0 to 1) vertically across the map, was used. For testing, this map had a size of 101 × 101 pixels. This ensures all combinations of oxygen saturation and total hemoglobin abundance are evenly represented in the map (note, this is a larger version of the phantom used in optimization, which was 21 × 21 pixels) [ Fig. 3(a)]. Test images with various degrees of gaussian white noise added to the signal hypercube and various degrees of gaussian white noise added to the simulated images were also generated to assess the effect of noise on filter set performance.
For qualitative assessment of the performance, a ground truth map inspired by biological images of vascular networks was prepared [ Fig. 3(b)]. This included a radial network of vessels with decreasing vessel diameter at the periphery of the image. Inside these vessels the total hemoglobin abundance was set to a value of 1. Outside vessels, the total hemoglobin abundance was set to a value of 0.5 to represent regions where the signal contribution from small unresolvable capillaries is mixed with the non-vascular background, ultimately representing a homogenous tissue background.

Results
Spectral filter sets were chosen using four methods: evenly spaced filters, like those found in an 'off-the-shelf' system (1); filters optimized by minimizing mutual information between filters (2); and filters optimized by minimizing the RMSE in hemoglobin abundance prediction, using gradient descent (3); and using a genetic algorithm (4). The selected filters for n = 2, 3, 4, 5, 9, 16 and 25 are shown in Fig. 4 for mutual information, gradient descent, and genetic algorithms. The time taken for each method is shown in Table S1. The selected filters for all n are shown in Fig. S2. The hierarchical clustering of filters selected by minimizing mutual information is evident in Fig. 4(a). A heat map of mutual information is shown in Fig. S3. The filters optimized by gradient descent [Fig. 4(a)] have center wavelengths close to the evenly spaced starting center wavelengths, suggesting the algorithm became trapped in local minima. In contrast, the filters optimized by genetic algorithm [Fig. 4(c)] are more unevenly distributed suggesting a better global optimization. The performance of each filter set was determined by simulating imaging of an independent test hypercube generated using the gradient abundance map shown in Fig. 3(a). Following spectral unmixing, the root-mean-square-error (RMSE) in estimated abundance [ Fig. 5(a)] and SO 2 [ Fig. 5(b)] were calculated. For high numbers of filters, n > 15, there is little absolute difference in performance between filter sets selected using each of the four methods, as most of the wavelength range is sampled. However, for sparse sampling (small numbers of filters; n < 10), filter sets optimized using a genetic algorithm result in superior performance to filters optimized by all other approaches, with up to a 37% reduction in RMSE-SO 2 compared to evenly spaced filters [ Fig. 5(c)]. For example, optimization of n = 3 filters using the genetic algorithm results in a 39% decrease in RMSE-abundance (0.06 vs. 0.10) and a 29% decrease in RMSE-SO 2 (0.064 vs. 0.090) compared with evenly spaced filters.
Optimization based on the alternative merit functions (Section 2.7) that maximize spectral angle (SA) and root-mean-square-difference (RMSD) between HbO 2 and Hb spectra perform poorly at low n, but better at n > 7, as do filters optimized using mutual information, as can be seen in the percentage reduction of RMSE-SO 2 versus evenly spaced filters [ Fig. 5(c)], but at high filter numbers, this absolute improvement is small [ Fig. 5(b)].
The error in estimated SO 2 is shown in Fig. 6 for evenly spaced filter sets and optimized filter sets for n = 3, 4, 9 and 25. In all cases the error in estimated SO 2 is largest at low THb abundance due to lower signal to noise ratio. For n = 3 and n = 4, the error in SO 2 for intermediate A THb ∼0.3-0.6 is significantly lower using optimized filter sets than using evenly spaced filter sets.
To assess the effects of noise, gradient test images [ Fig. 3(a)] with various degrees of gaussian white noise added to the signal hypercube and various degrees of gaussian white noise added to the simulated images were generated. These images were 'imaged' with the optimized filter sets for n = 3, unmixed, and root-mean-square-error in estimated SO 2 (RMSE-SO 2 ) calculated. The results are shown in Fig. 7. Filter sets optimized to minimize RMSE via genetic algorithm remained the highest performing filter sets across noise conditions. Figure 8 shows the estimated abundance maps for imaging with n = 3, 4, 9 and 25 filter sets selected using three different methods: evenly spaced filters, filters minimizing mutual information and filters optimized using a genetic algorithm (GA). Further examples are shown in Fig. S4. For n = 16 and n = 25, images are similarly accurate for filter sets chosen by all methods. For n = 3 and n = 4, the images simulated using optimized filters are visibly less noisy and more accurate than the images simulated using evenly spaced filters and filters selected to minimize 5 10 15 20 25  The ground truth abundance maps for two regions of the vessel image are shown (top). For each of these regions, the estimated abundance maps are shown for imaging with n = 3, 4, 9 and 25 filters selected using three different methods: evenly spaced filters, filters optimized by minimizing mutual information and filters optimized by minimizing root-mean-squareerror (RMSE) via a genetic algorithm (GA). The red channel represents the abundance of oxyhemoglobin. The blue channel represents the abundance of deoxyhemoglobin.

Discussion
The results clearly demonstrate the merit of tailoring spectral filter sets towards specific biomedical problems. The RMSE in determining SO 2 was 24-37% lower when imaging with n = 2-10 optimized spectral filters compared with the same numbers of evenly spaced filters. Moreover, the results make clear the importance of performing this optimization with the end-goal clinical challenge in mind by using appropriate functions of merit to assess filter set performance, in this case the RMSE in estimating hemoglobin abundances. Filter sets optimized using a merit function based on accurate estimation of hemoglobin abundance resulted in significantly more accurate estimation of SO 2 compared with spectral filter sets chosen to minimize mutual information, maximize spectral angle between HbO 2 and Hb spectra or maximize root-mean-square-difference between HbO 2 and Hb spectra (69-79%, 43-80% and 43-67% decrease in RMSE-SO 2 for imaging with n = 2-5 filters chosen to minimize mutual information, maximize spectral angle between HbO 2 and Hb spectra or maximize root-mean-square-difference between HbO 2 and Hb spectra respectively).
The results also show that optimized filter sets consistently perform better than non-optimized or poorly optimized filter sets across a wide range of noise levels. The noise levels in experimental images will depend on the experimental setting as well as the acquisition parameters such as exposure time and illumination power. In standard practice, we can assume these would be adjusted by the user to achieve adequately small noise levels and ensure an optimized filter set's performance matches that seen in simulation. Otherwise, we suggest users perform their own optimization using noise levels tuned to those expected in their own experimental setups. Now is the time to perform spectral filter optimizations. Many spectral imaging approaches are on the cusp of clinical translation, but despite the number of publications, very few in vivo studies have been conducted [27]. Commercial devices for hyperspectral imaging of oxygenation, such as the TIVITA series imagers (Diaspective Vision GmbH, Germany), HyperView (Hypermed Imaging, Inc., USA) and Snapshot NIR (Kent Imaging, Inc., Canada), are rapidly emerging, and multiple hyperspectral imaging start-ups have been founded in the past few years. It is our belief that wider adoption of proper spectral band selection will result in more effective multispectral imaging systems and consequently, more widespread adoption of multispectral imaging technologies in clinic.
We foresee a future where ubiquitous 'one-size-fits-all' RGB imaging systems in endoscopy, laparoscopy and surgical microscopy are replaced with application-specific multispectral imaging systems tailored to the measurement of biological signatures relevant to the said application. Proponents of the 'one-size-fits-all' approach might point out the increased cost associated with purchasing and maintaining multiple application-specific imaging tools, but this is already the norm in surgical instruments, where a plethora of application-specific scalpels, forceps, scissors, retractors, and clamps have long been available. Still, the relatively high cost of imaging sensors might previously have disqualified this comparison, but recent innovations in optical filter fabrication [42,43] increasingly allow for low-cost manufacture of custom sensors.
An alternative and promising approach to multispectral imaging is the use of sequential spectral filtering of the illumination, allowing the use of common monochrome detectors, as is already in use in flexible endoscopy for RGB, autofluorescence imaging and NBI. The increasing availability of high-power narrow-band LEDs might facilitate the use of multi-LED arrays for illumination. This would enable the continued use of ubiquitous 'one-size-fits-all' hardware, with software instructing the multi-LED array to provide application-specific sequential illumination as determined by spectral filter optimization per indication.
We have identified several limitations to this work. To ensure applicability across acquisition approaches, the present study did not consider the trade-off between increasing the number of filters and decreasing spatial and/or temporal resolution, but this has been explored previously for multispectral filter arrays. Ultimately, the trade-off between spectral and temporal resolution will be a somewhat subjective choice that depends on the intended end-users and application. Another limitation of the current study is the use of standard fitting by minimization of a sum of square errors cost function to fit the empirical model to the simulated images and determine estimated abundances. Increasingly, convolutional neural networks (CNNs) are being used to improve the interpretation of biomedical images. An embedded CNN-based approach to spectral band selection should be explored for use in biomedical spectral band optimization once appropriate labelled datasets are available [44]. Finally, the present study used a simple empirical model of reflectance; future work may expand this to a 3D tissue model using Monte Carlo photon transport simulations. Nevertheless, the empirical model's good fit to experimental esophageal tissue data is encouraging. For future work, we advise users to employ experimentally acquired ground truth hypercubes as inputs, so optimization is based on the most accurate spectral properties of the tissue of interest.
In the present study we chose to focus on Gaussian filter profiles as these can be used to approximate the spectral filter profiles of narrow-band light sources (e.g. LED arrays), tunable filters (e.g. LCTFs, AOTFs) and spectrally resolved detector arrays (e.g. micro-pixelated filters). Still, other users may require optimization of different filter shapes (e.g. dichroic band pass). The presented algorithm is designed such that alternative filter shapes could be added as inputs to the optimization, and we encourage users to perform optimizations to suit their needs. The results of the present study demonstrate the potential of spectral band optimization, and we hope this will persuade researchers to use the presented algorithms, with their own inputs, to enable optimization tailored towards their own specific applications.

Conclusions
We demonstrate the optimization of spectral filter sets enables more accurate estimation of oxygen saturation in spectral imaging, with up to a 37% reduction in root-mean-square-error compared with the use of generalist sensors. This work clearly demonstrates the merit of tailoring spectral filter sets towards specific biomedical problems. Wider adoption of this approach will result in more effective multispectral image systems and consequently, more widespread adoption of multispectral imaging technologies in clinic.