1. Introduction {#sec1}
===============

Multispectral reflectance is a broadly established method for estimating concentrations of multiple absorbers in tissue. It is often used to decouple and estimate quantities of oxyhemoglobin (${HbO}_{2}$), deoxyhemoglobin (HbR), and water ($H_{2}O$), among other absorbers, through spectral fitting based on least squares, provided the acquisition of multiple wavelengths is available (a minimum of three wavelengths is needed for ${HbO}_{2}$, HbR, and $H_{2}O$). The wavelengths acquired have a considerable effect on the accuracy of the concentration estimates that are generated. In the approach presented here, wavelengths are selected by maximizing the product of the singular values of the scattering-modulated absorption matrix. Analysis shows that this selection method is robust and stable over a wide range of optical properties and that it improves the accuracy of the resulting estimates relative to algorithms previously published. Finally, a practical application where optimal light-emitting diodes (LEDs) are chosen to estimate blood concentrations is described.

2. Methods {#sec2}
==========

2.1. Spectral Fitting {#sec2.1}
---------------------

Concentrations of multiple absorbers in tissue can be estimated from its wavelength-dependent reflectance signal. For each wavelength, the fraction of reflectance collected can be modeled with Beer's law $$R(\lambda) = \exp\lbrack - \mu_{a}(\lambda)L(\lambda)\rbrack,$$where $\mu_{a}(\lambda)$ is the wavelength-dependent absorption coefficient represented in the medium as a linear combination of the absorption coefficients of all components present. To illustrate the benefits of our approach, we consider oxyhemoglobin, deoxyhemoglobin, and water, and, accordingly, the absorption coefficient is $$\mu_{a}(\lambda) = BS\mu_{oxy}(\lambda) + B(1 - S)\mu_{deoxy}(\lambda) + W\mu_{\text{water}}(\lambda),$$where $L(\lambda)$ is the average path length of photons in the medium at wavelength $\lambda$, which varies according to the absorption ($\mu_{a}$) and scattering ($\mu_{s}^{\prime}$) coefficients of the tissue, and is computed using an adaption of the diffusion model described by Farrell et al.[@r1] $B$, $S$, and $W$ quantify the blood concentration, oxygen saturation, and water concentration, respectively. Parameter estimation can be refined through more advanced models. For example, adding specific layers of tissue can be advantageous if the laminar distribution of chromophores is known (e.g., in the case of epidermial melanin). Additionally, correction factors accounting for vessel size can be added depending on the underlying structures that are being considered (notably, blue light does not penetrate to the center of large vessels and, therefore, underestimates blood content). The optical depth is defined as the natural logarithm of the reflectance signal and is needed to extract information about concentration through spectral fitting $$- \ln\, R(\lambda) = \mu_{a}(\lambda)L(\lambda).$$

Given an optical depth equation for each wavelength, the problem is decomposed into a matrix system $$\frac{\underbrace{\begin{bmatrix}
{- \ln\lbrack R(\lambda_{1})\rbrack} \\
{- \ln\lbrack R(\lambda_{2})\rbrack} \\
 \vdots \\
{- \ln\lbrack R(\lambda_{N})\rbrack} \\
\end{bmatrix}}}{\mathbf{OD}} = \frac{\underbrace{\begin{bmatrix}
{\mu_{a,{oxy}}(\lambda_{1})L(\lambda_{1})} & {\mu_{a,{deoxy}}(\lambda_{1})L(\lambda_{1})} & {\mu_{a,\text{water}}(\lambda_{1})L(\lambda_{1})} \\
{\mu_{a,{oxy}}(\lambda_{2})L(\lambda_{2})} & {\mu_{a,{deoxy}}(\lambda_{2})L(\lambda_{2})} & {\mu_{a,\text{water}}(\lambda_{2})L(\lambda_{2})} \\
 \vdots & \vdots & \vdots \\
{\mu_{a,{oxy}}(\lambda_{N})L(\lambda_{N})} & {\mu_{a,{deoxy}}(\lambda_{N})L(\lambda_{N})} & {\mu_{a,\text{water}}(\lambda_{N})L(\lambda_{N})} \\
\end{bmatrix}}}{\mu_{a}L}\frac{\underbrace{\begin{bmatrix}
{BS} \\
{B(1 - S)} \\
W \\
\end{bmatrix}}}{\mathbf{C}},$$where $\mathbf{OD}$ is a $1 \times N$ vector of optical depths, $\mu_{a}L$ is a $N \times 3$ matrix of pathlength-modulated absorption coefficients associated with the absorbers of interest, and $\mathbf{C}$ is the concentration vector. The concentration vector can be computed using the Moore--Penrose pseudoinverse.[@r2] Positivity constraints are added to ensure that no negative concentrations will be computed,[@r3] yielding the least square fit $$\mathbf{C} = \frac{\underbrace{{\lbrack{(\mu_{a}L)}^{T}(\mu_{a}L)\rbrack}^{- 1}{\lbrack\mu_{a}L\rbrack}^{T}}\mathbf{OD}}{\text{Moore} - \text{Penrose\ pseudoinverse}}.$$

In principle, any three wavelengths can be used to generate concentration estimates, but their selection affects the resulting accuracy dramatically.[@r4][@r5][@r6][@r7][@r8][@r9]^--^[@r10]

2.2. Wavelength Selection {#sec2.2}
-------------------------

In this paper, optimal wavelengths are selected from an oversampled set of wavelengths and iteratively removed until the desired (minimum) number of wavelengths is found, based on a selection algorithm used by Luke et al.[@r4]^,^[@r5] For each row (i.e., wavelength) being removed, selection criteria are computed to evaluate the importance of that wavelength. For clarity, a pseudocode for the selection algorithm is outlined here:

1.  Start with matrix $\mu_{a}L$ that has an oversampled number of wavelengths.

2.  Remove the first row (i.e., wavelength of matrix $\mu_{a}L$). Compute the value of the selection metric for this reduced matrix.

3.  Iteratively repeat step 2 for all rows of $\mu_{a}L$.

4.  Eliminate the row whose removal yielded the largest/smallest selection metric (depending on the metric).

5.  Repeat these steps until the desired number of wavelengths is achieved.

Selection of optimal wavelengths is achieved by maximizing or minimizing a selection metric. This metric is often based on the singular values of the system's absorption matrix (in this case, $\mu_{a}L$). These singular values indicate the orthogonality of the spectra being fitted and can be found through singular value decomposition. Previous studies have investigated the selection of optimal wavelengths for a wide variety of spectroscopic applications, including photoacoustic imaging,[@r4]^,^[@r5] spatial frequency-domain imaging,[@r6] and optical tomography.[@r8]^,^[@r9] Most of these reports base selection on the matrix condition number (ratio of the largest to the smallest singular value).[@r4]^,^[@r6]^,^[@r8] Luke and Emelianov[@r5] also maximized the smallest singular value, which prevents loss in rank, and reported improved results using this approach. In this paper, we maximized the product of all singular values, which yielded concentration estimates that were more accurate than those generated from wavelengths selected by previously published algorithms. While the condition number indicates the spread in singular values and using the smallest singular values prevents the system from losing rank, both approaches discard information contained in the rest of the singular values whereas the product of singular values senses information from the full matrix spectrum.

To evaluate the performance of sets of optimal wavelengths generated by various methods, simulated spectra were generated from randomly chosen optical properties and concentrations of ${HbO}_{2}$, HbR, and $H_{2}O$ from ranges reported in [Table 1](#t001){ref-type="table"}. Ten percent (10%) white Gaussian noise was added to each of these simulated spectra, and least square fits were performed to estimate concentrations. Root mean squared (RMS) errors between the estimates and the true concentrations were recorded. This process was repeated $10^{3}$ times for 10 randomly selected optical properties, and the resulting averaged RMS error was used to compare algorithms and evaluate the quality of selected wavelengths.

###### 

Ranges of the parameters used to generate random absorption spectra.

  Parameter   $B$              $S$           $W$               $a$     $b$
  ----------- ---------------- ------------- ------------- ----------- --------
  Range       0.001 to 0.003   0.5 to 0.95   0.65 to 0.8    0.33 to 3  0 to 3

Note: $B$ is the blood volume fraction, $S$ is the average arteriovenous oxygen saturation of hemoglobin, $W$ is water volume fraction, and $a$ and $b$ specify the reduced scattering coefficient $\mu_{s}^{\prime} = a{(\lambda/500\text{  }{nm})}^{- b}$, such that $a = \mu_{s}^{\prime}(500\text{  }{nm})$

3. Results {#sec3}
==========

3.1. Selection Criterion {#sec3.1}
------------------------

To demonstrate that the center singular value of $\mu_{a}L$ has an important impact on the quality of the resulting least square estimates, simulated spectra of random optical properties were generated to populate $\mu_{a}L$ matrices, and their center singular values were artificially increased or decreased within the range of their extreme singular values. Least square fitting was performed with the scaled center singular value, and RMS errors were computed and are shown in [Fig. 1](#f1){ref-type="fig"}.

![Average relative RMS error generated when performing spectral fitting for ${HbO}_{2}$, HbR, and $H_{2}O$ with a center singular value varying from the smallest to largest singular value of the system.](JBO-023-071202-g001){#f1}

The error decrease shows that the size of the center singular value has a considerable effect on the accuracy of the concentration estimates. Artificially modifying the center singular value maintained both the condition number and smallest singular value of the system as constants whereas the average RMS errors varied over orders of magnitude. Accordingly, maximizing the product of singular values when selecting optimal wavelengths is a simple but desirable selection criterion. To further illustrate the advantage of this selection metric, average RMS errors generated when fitting with wavelengths chosen by different selection criteria were compared in [Fig. 2](#f2){ref-type="fig"}. As a control, the RMS error associated with wavelengths selected randomly and linearly within the 480- to 1000-nm range was also computed.

![Average relative RMS error of concentration estimates generated by fitting ${HbO}_{2}$, HbR, and $H_{2}O$ for wavelengths chosen using different selection methods. RMS error was averaged over 100 spectra of random optical properties contained in [Table 1](#t001){ref-type="table"} ranges. The methods compared are: (1) minimizing the condition number of $\mu_{a}L$, (2) maximizing the product of the singular values, (3) maximizing the minimum singular value,[@r5] (4) a random choice of wavelengths, and (5) linear spacing of wavelengths.](JBO-023-071202-g002){#f2}

These results demonstrate that using the product of the singular values yields improved results over condition number or smallest singular value, especially when selecting a larger number of wavelengths. Interestingly, random and linearly spaced wavelength selection outperformed those based on condition number and the smallest singular value. Methods that rely on the smallest singular values are affected by spectra that are highly disproportioned (as is the case here, with water having an absorption coefficient that is orders of magnitude lower than blood).

3.2. Optimal Wavelengths {#sec3.2}
------------------------

Using our selection algorithm based on maximizing the product of singular values, the three optimal wavelengths for estimating blood concentration and saturation in tissue were computed for a wide variety of optical properties. For the majority (76%) of optical property combinations tested, the same trio of optimal wavelengths was selected. These wavelengths are shown in [Fig. 3](#f3){ref-type="fig"}.

![Optimal wavelengths selected by maximizing the product of the singular values of the system.](JBO-023-071202-g003){#f3}

In the other cases, which were mostly when oxygen saturation was lower than 60%, the center wavelength selected switched from 604 to 596 nm. Variations in the wavelengths picked are expected when modifying optical properties since they alter the shape of the absorption spectra used to compute the optimal wavelengths. Even if, for some combinations of optical properties, the selection algorithm found wavelengths that differed slightly from the ones in [Fig. 3](#f3){ref-type="fig"}, these choices generated similar RMS errors. The concentration estimates generated with the set of optimal wavelengths highlighted in [Fig. 3](#f3){ref-type="fig"} yielded RMS errors that were smaller than the ones generated with other wavelength trios, on average. The region of low RMS error was constant for a radius of about 10 nm around the optimal wavelengths, as shown in [Fig. 4](#f4){ref-type="fig"}, in which the change in RMS error relative to the error generated with the optimal wavelengths is presented.

![Relative change in average RMS error for wavelength trios around the optimal wavelengths selected by maximizing the product of the singular values. Each square represents a 5-nm variation.](JBO-023-071202-g004){#f4}

Interestingly, some wavelengths can be moved with greater freedom than others. For example, increasing the 578-nm wavelength changes RMS error nominally, while decreasing it produces RMS errors that are orders of magnitude greater.

3.3. Additional Chromophores {#sec3.3}
----------------------------

Depending on the tissue being imaged, more (or different) absorbers may need to be considered. They can easily be included in the selection process by adding their associated spectra to the matrix of pathlength-modulated bases ($\mu_{a}L$) used to compute the singular values. As a demonstration, the spectra of ICG (for a concentration of $1290\text{  }\mu M$), lipids, and a linear offset (which removes linear background sources) were added to the basis spectra matrix, and 6 wavelengths were selected with these additional components. As in the 3 wavelength selection case, optical properties and concentrations were chosen randomly 10 times within the ranges outlined in [Table 1](#t001){ref-type="table"}, for concentrations of ICG and lipids varying between 0.02 and 0.06 (%). The optimal wavelengths found by our algorithm are presented in [Fig. 5](#f5){ref-type="fig"}. The average RMS error of the concentration estimates generated with these wavelengths was compared with the average errors generated with other algorithms, and the results are summarized in [Table 2](#t002){ref-type="table"}.

![Wavelengths chosen for decoupling ${HbO}_{2}$, HbR, $H_{2}O$, ICG, lipids, and a linear offset.](JBO-023-071202-g005){#f5}

###### 

Average RMS error for fitting six components.

  Selection method    Linear spacing   Random selection   Condition number   Smallest singular value   Product of singular value
  ------------------- ---------------- ------------------ ------------------ ------------------------- ---------------------------
  Average RMS error   0.39             1.45               0.37               0.30                      0.13

Once again, selecting wavelengths based on the product of singular values yielded the smallest RMS error and outperformed previously published algorithms.

4. Practical Application {#sec4}
========================

The selection algorithm outlined in this paper is easily adapted to multiple practical situations. As an example, it is used to select the 3 best LEDs for fitting ${HbO}_{2}$, HbR, and $H_{2}O$ from an inventory of 20 possibilities. The input set of LEDs was taken from Gareau.[@r11] Although LED emission spectra usually have varying amplitudes, they have been normalized for this analysis to account only for their emission wavelength. The normalized spectra are shown in [Fig. 6](#f6){ref-type="fig"}.

![Emission spectra of 20 LEDs considered for diffuse reflectance acquisition of ${HbO}_{2}$, HbR, and $H_{2}O$.](JBO-023-071202-g006){#f6}

LEDs were selected by the algorithm outlined in this paper, by replacing the $\mu_{a}L$ matrix with the following LED matrix $$\begin{bmatrix}
{\int\mu_{a,{oxy}}L\,{LED}_{1}d\lambda} & {\int\mu_{a,{deoxy}}L\,{LED}_{1}d\lambda} & {\int\mu_{a,{water}}L\,{LED}_{1}d\lambda} \\
{\int\mu_{a,{oxy}}L \cdot {LED}_{2}d\lambda} & {\int\mu_{a,{deoxy}}L\,{LED}_{2}d\lambda} & {\int\mu_{a,{water}}L\,{LED}_{2}d\lambda} \\
 \vdots & \vdots & \vdots \\
{\int\mu_{a,{oxy}}L\,{LED}_{20}d\lambda} & {\int\mu_{a,{deoxy}}L\,{LED}_{20}d\lambda} & {\int\mu_{a,{water}}L\,{LED}_{20}d\lambda} \\
\end{bmatrix},$$where ${LED}_{i}$ is the emission spectrum of LED $i$. This matrix considers the bandwidth of the LED's emission as well as the modulated absorption spectra of the components fitted.

Executing the selection algorithm on this matrix returned the 3 best LEDs: LED\# 9, LED\# 17, and LED\# 20 centered at 496, 668, and 935 nm, respectively, and yielded wavelength peaks having slight variations from the optimal wavelengths picked by our algorithm. This outcome is explained by the fact that, when dealing with LEDs, not all wavelengths are available and that other factors, such as LED bandwidth, are considered.

5. Conclusion {#sec5}
=============

This paper presented and evaluated a wavelength selection criterion that produced improved results relative to previously published strategies by maximizing the product of singular values of the scattering-modulated absorption matrix. Using this approach, three optimal wavelengths were found for diffuse reflectance spectral fitting of ${HbO}_{2}$, HbR, and $H_{2}O$. These optimal wavelengths proved to be stable over a wide range of optical absorption and scattering properties. The method is easily extended to more and different absorbers and achieved similarly improved results. Modifications to the algorithm were demonstrated to account for LED broadband emission spectra instead of distinct discrete wavelengths, and the process was used to select LEDs that are optimal for fitting blood components. While this paper prioritized the diffuse reflectance fitting of blood, the algorithm is readily applied to applications that require spectral fitting of other absorbers and/or fluorophores.

This work was supported in part by the National Institutes of Health under Grant No. R01NS052274 The authors would like to thank Geoffrey P. Luke, Jonathan D. Olson, and Jaime B. Bravo for helpful discussions on the subject of wavelength selection.

The authors have no conflicts of interests to disclose.

**Mikael Marois** received his BS degree in engineering physics from the Polytechnique Montreal in 2015. Currently, he is a graduate student working with the Optics in Medicine group, the Thayer School of Engineering, Dartmouth College, Hanover, New Hampshire, USA. His research is focused on interpreting optical signals measured during intraoperative neurosurgeries.

**Steven L. Jacques** received his BS degree in biology from MIT in 1972, his MSEE degree in electrical engineering and computer science, and his PhD in biophysics and medical physics from the University of California, Berkeley in 1984. He has served at Massachusetts General Hospital/Harvard Medical School, the University of Texas M.D. Anderson Cancer Center, the Oregon Health and Science University, and currently at Tufts University. His research focuses on light propagation and interactions with tissues.

**Keith D. Paulsen** is a professor of biomedical engineering at Thayer School of Engineering and of radiology and surgery at Dartmouth--Hitchcock Medical Center. He is an expert in biomedical imaging and computational modeling. His research has focused on the development and translation of advanced imaging technology, primarily for cancer detection, diagnosis, therapy monitoring, and surgical guidance. He has authored more than 350 publications and has maintained an active research program, continuously funded by the National Institutes of Health for over 25 years.
